NYC COVID-19 Social Disparities
Carrión et al.
2020-11-24
- Session Configuration
- Functions
- Load Data
- Census Data
- Grocery stores per area
- Subway station locations
- Residential area
- COVID Tests
- Create data frames of all above information
- Tract to ZCTA weighted assignment
- Part 1: Creation of BWQS Neighborhood Infection Risk Scores
- BWQS by Tract
- Part 2: Compare capacity to socially distance (as measured by transit data) by neighborhood-level risk scores
- Part 3: Spatial analysis of mortality in relation to BWQS scores
- Subway analyses
- Appendix
t_start = Sys.time()# CRAN packages:
library(tidyverse)
library(sf)
library(lubridate)
library(tidycensus)
library(ggExtra)
library(ggridges)
library(ggsn)
library(ragg)
library(rstan)
library(drc)
library(spdep)
library(broom)
library(MASS)
library(spatialreg)
library(here)
library(pdftools)
library(matrixStats)
library(egg)
library(ggpubr)
library(scales)
library(qs)
library(corrplot)
library(readxl)
library(splines)
library(magic)
library(httr)
library(jsonlite)
library(DHARMa)
library(kableExtra)
library(lwgeom)
# Github packages available via:
# remotes::install_github("justlab/Just_universal")
# remotes::install_github("justlab/MTA_turnstile")
library(Just.universal) Session Configuration
#### SESSION CONFIGURATION ####
options(dplyr.summarise.inform=FALSE)
# data will default to a subfolder "data/" within working directory
# unless 1. set by an environment variable:
data.root = Sys.getenv("COVID_DATA")
# or 2. set with an alternative path here:
if (data.root == "") data.root = here("data")
if (!dir.exists(data.root)) dir.create(data.root)
print(paste("data being downloaded into directory", dQuote(data.root)))## [1] "data being downloaded into directory \"/data/scratch/COVID_runtime/data\""# Get some data from the git repository rather than downloading from original
# source, to avoid changes in model results due to updated data
use_repo_data = TRUE
if(Sys.getenv("MTA_TURNSTILE_DATA_DIR") == ""){
# set up default download location for MTA turnstile data
mta_dir = file.path(data.root, "mta_turnstile")
if(!dir.exists(mta_dir)) dir.create(mta_dir, recursive = TRUE)
Sys.setenv(MTA_TURNSTILE_DATA_DIR = mta_dir)
}
library(MTA.turnstile)## MTA.turnstile data directory: /data/scratch/COVID_runtime/data/mta_turnstilet_turnstile_1 = Sys.time()
# output path for figures
fig.path = here("figures")
if(!dir.exists(fig.path)) dir.create(fig.path)
export.figs = TRUE
if(export.figs) message("Saving figures to:\n ", fig.path) else message("Not saving figures")## Saving figures to:
## /data/scratch/COVID_runtime/figures# To generate census data, you need an API key, which you can request here: https://api.census.gov/data/key_signup.html
#census_api_key("INSERT YOUR CENSUS API KEY HERE", install = TRUE)
if(Sys.getenv("CENSUS_API_KEY")=="") warning("Census API Key Missing")
# pairmemo function cache
pairmemo.dir = file.path(data.root, "pairmemo")
dir.create(pairmemo.dir, showWarnings = F)
pm = function(...) pairmemo(
directory = pairmemo.dir,
n.frame = 2,
...)About pairmemo:
pairmemo caches the results of functions on disk based on their input parameters.
If you run this script and then make changes to code, we recommend that you clear the pairmemo cache of any functions that run in the lines following your change before running the script again.
The following commented lines would clear the pairmemo cache for each function. If you do not change the input Census or Pluto data, you would likely only need to run the last one.
#pairmemo.clear(get.Pluto)
#pairmemo.clear(get.qn.blocks)
#pairmemo.clear(acs.main)
#pairmemo.clear(get_essential_acs)
#pairmemo.clear(get.tract.res)
#pairmemo.clear(fit_BWQS_model)Functions
#### Functions ####
read_w_filenames <- function(flnm) {
read_csv(flnm) %>%
mutate(filename = flnm)
}
extract_waic <- function (stanfit){
log_lik <- rstan::extract(stanfit, "log_lik")$log_lik
dim(log_lik) <- if (length(dim(log_lik)) == 1)
c(length(log_lik), 1)
else c(dim(log_lik)[1], prod(dim(log_lik)[2:length(dim(log_lik))]))
S <- nrow(log_lik)
n <- ncol(log_lik)
lpd <- log(colMeans(exp(log_lik)))
p_waic <- colVars(log_lik)
elpd_waic <- lpd - p_waic
waic <- -2 * elpd_waic
loo_weights_raw <- 1/exp(log_lik - max(log_lik))
loo_weights_normalized <- loo_weights_raw/matrix(colMeans(loo_weights_raw),
nrow = S, ncol = n, byrow = TRUE)
loo_weights_regularized <- pmin(loo_weights_normalized, sqrt(S))
elpd_loo <- log(colMeans(exp(log_lik) * loo_weights_regularized)/colMeans(loo_weights_regularized))
p_loo <- lpd - elpd_loo
pointwise <- cbind(waic, lpd, p_waic, elpd_waic, p_loo, elpd_loo)
total <- colSums(pointwise)
se <- sqrt(n * colVars(pointwise))
return(list(waic = total["waic"], elpd_waic = total["elpd_waic"],
p_waic = total["p_waic"], elpd_loo = total["elpd_loo"],
p_loo = total["p_loo"]))
}
# Download a file, update metadata records, and load it with function `f`
# File metadata is stored in a sqlite file, by default in data/downloads/meta.sqlite
download = function(url, to, f, ...){
f(download.update.meta(url, file.path(data.root, "downloads"), to),
...)
}Load Data
#### Load Data ####Progress bars are a single line when run interactively, but print every refresh with knitr
# Subway ridership data
Subway_ridership_by_UHF <- relative.subway.usage(2020L, "nhood")## Reading turnstile files## Processing## Parsing timestamps## Sorting## Setting station names## Dropping 3 unlocated stations## Dropping 278,230 observations## Decumulating - entries##
|
| | 0%
|
| | 1%
|
|= | 1%
|
[manually edited for length]
|
|======================================================================| 99%
|
|======================================================================| 100%## Dropped 189 sources of 5,292## Dropped 1,102,556 observations of 60,962,516## Decumulating - exits##
|
| | 0%
|
| | 1%
|
|= | 1%
|
[manually edited for length]
|
|======================================================================| 99%
|
|======================================================================| 100%## Dropped 173 sources of 5,292## Dropped 868,952 observations of 60,962,516t_turnstile_2 = Sys.time()
# get the Pluto dataset from #https://www1.nyc.gov/site/planning/data-maps/open-data/dwn-pluto-mappluto.page
pm(fst = T,
get.Pluto <- function() download(
"https://www1.nyc.gov/assets/planning/download/zip/data-maps/open-data/nyc_pluto_20v3_csv.zip",
"pluto.zip",
function(p)
read_csv(unz(p, "pluto_20v3.csv"), col_types = cols(spdist2 = col_character(),
overlay2 = col_character(),
zonedist4 = col_character()))[,
c("landuse", "bbl", "numfloors", "unitstotal", "unitsres",
"zipcode")]))
Pluto <- as.data.frame(get.Pluto())
# If you get an error in running this function, the download may be incomplete/corrupt
# In that case, delete the pluto.zip file and try again:
#unlink(file.path(data.root, "downloads", "pluto.zip"))
if (file.exists(file.path(data.root, "Bldg_Footprints.qs"))) {
Bldg_Footprints <- qread(file.path(data.root, "Bldg_Footprints.qs"))
} else {
Bldg_Footprints <- download(
# https://data.cityofnewyork.us/Housing-Development/Building-Footprints/nqwf-w8eh
"https://data.cityofnewyork.us/api/geospatial/nqwf-w8eh?method=export&format=Shapefile",
"building_footprints.zip",
function(p)
st_read(paste0("/vsizip/", p)))
qsave(Bldg_Footprints, file.path(data.root, "Bldg_Footprints.qs"))
}## Multiple layers are present in data source /vsizip//data/scratch/COVID_runtime/data/downloads/building_footprints.zip, reading layer `geo_export_aad3a557-b481-4681-b754-34463d92a741'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.## Warning in evalq((function (..., call. = TRUE, immediate. = FALSE, noBreaks. =
## FALSE, : automatically selected the first layer in a data source containing more
## than one.## Reading layer `geo_export_aad3a557-b481-4681-b754-34463d92a741' from data source `/vsizip//data/scratch/COVID_runtime/data/downloads/building_footprints.zip' using driver `ESRI Shapefile'
## Simple feature collection with 1084927 features and 15 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -74.2555 ymin: 40.49843 xmax: -73.70006 ymax: 40.91505
## CRS: 4326ZCTA_by_boro <- download(
"https://www.health.ny.gov/statistics/cancer/registry/appendix/neighborhoods.htm",
"uhf_neighborhoods.html",
function(p)
{# XML::readHTMLTable doesn't identify the columns correctly.
x = str_match_all(read_file(p), regex(dotall = T, paste0(
'<td headers="header1"[^>]+>\\s*(.+?)</td>',
'(.+?)',
'(?=<td headers="header1"|</table>)')))[[1]]
do.call(rbind, lapply(1 : nrow(x), function(i)
data.frame(boro = x[i, 2], zip = as.integer(
str_extract_all(x[i, 3], "\\b\\d{5}\\b")[[1]]))))})
# Download the specific day of test results by ZCTA being used
ZCTA_test_download <- download(
"https://raw.githubusercontent.com/nychealth/coronavirus-data/6d7c4a94d6472a9ffc061166d099a4e5d89cd3e3/tests-by-zcta.csv",
"2020-05-07_tests-by-zcta.csv",
identity
)
# Download COVID-19 testing data
ZCTA_test_series <- ZCTA_test_download %>%
map_df(~read_w_filenames(.)) %>%
mutate(date = as.Date(str_extract(filename, "[:digit:]{4}-[:digit:]{2}-[:digit:]{2}"))) %>%
dplyr::select(-filename)## Parsed with column specification:
## cols(
## MODZCTA = col_double(),
## Positive = col_double(),
## Total = col_double(),
## zcta_cum.perc_pos = col_double()
## )# UHF definitions by zip codes
UHF_ZipCodes <- UHF_ZipCodes <- download(
"http://www.infoshare.org/misc/UHF.pdf",
"uhf_zips.pdf",
function(p)
{x = str_match_all(pdf_text(p)[2],
"(\\d+)\\s+(\\S.+?\\S)\\s*([0-9,]+)")[[1]]
do.call(rbind, lapply(1 : nrow(x), function(i)
data.frame(code = x[i, 2], name = x[i, 3], zip = as.integer(
str_extract_all(x[i, 4], "\\b\\d{5}\\b")[[1]]))))})
# UHF shapefile
UHF_shp <- download(
"https://www1.nyc.gov/assets/doh/downloads/zip/uhf42_dohmh_2009.zip",
"nyc_uhf_nhoods_shapefile.zip",
function(p) read_sf(paste0("/vsizip/", p, "/UHF_42_DOHMH_2009")))
# NYC boroughs from NYC Open Data
NYC_basemap_shp <- download(
"https://data.cityofnewyork.us/api/geospatial/tqmj-j8zm?method=export&format=Shapefile",
"Borough_Boundaries.zip",
function(p){
unzip(p, exdir = file.path(data.root, "downloads"))
# open data platform generates a random UUID for every download
ufile = list.files(file.path(data.root, "downloads"), pattern = "geo_export_.*\\.shp", full.names = TRUE)
st_read(ufile, stringsAsFactors = FALSE, quiet = TRUE) %>% st_transform(., crs = 2263)
}
)
# DOHMH MODZCTA Shapefile
MODZCTA_NYC_shp <- download(
"https://data.cityofnewyork.us/api/geospatial/pri4-ifjk?method=export&format=Shapefile",
"Modified Zip Code Tabulation Areas (MODZCTA).zip",
function(p) read_sf(paste0("/vsizip/", p))
)
# Food outlets
if(use_repo_data){
food_retail <- read_csv(here("data", "retail_food_stores_2019-06-13.csv"))
} else {
food_retail <- download(
"https://data.ny.gov/api/views/9a8c-vfzj/rows.csv",
"retail_food_stores.csv",
read_csv)
}## Parsed with column specification:
## cols(
## County = col_character(),
## `License Number` = col_character(),
## `Operation Type` = col_character(),
## `Establishment Type` = col_character(),
## `Entity Name` = col_character(),
## `DBA Name` = col_character(),
## `Street Number` = col_character(),
## `Street Name` = col_character(),
## `Address Line 2` = col_logical(),
## `Address Line 3` = col_logical(),
## City = col_character(),
## State = col_character(),
## `Zip Code` = col_double(),
## `Square Footage` = col_double(),
## Location = col_character()
## )# Download deaths by ZCTA as of May 23rd
deaths_by23May2020_by_zcta <- download(
"https://raw.githubusercontent.com/nychealth/coronavirus-data/8d88b2c06cf6b65676d58b28979731faa10c193c/data-by-modzcta.csv",
"2020-05-23_data-by-modzcta.csv",
read_csv
)## Parsed with column specification:
## cols(
## MODIFIED_ZCTA = col_double(),
## NEIGHBORHOOD_NAME = col_character(),
## BOROUGH_GROUP = col_character(),
## COVID_CASE_COUNT = col_double(),
## COVID_CASE_RATE = col_double(),
## POP_DENOMINATOR = col_double(),
## COVID_DEATH_COUNT = col_double(),
## COVID_DEATH_RATE = col_double(),
## PERCENT_POSITIVE = col_double()
## )# download MODZCTA to ZCTA crosswalk, current version from repo
modzcta_to_zcta <- download(
"https://raw.githubusercontent.com/nychealth/coronavirus-data/master/Geography-resources/ZCTA-to-MODZCTA.csv",
"ZCTA-to-MODZCTA.csv",
read_csv
)## Parsed with column specification:
## cols(
## ZCTA = col_double(),
## MODZCTA = col_double()
## )modzcta_to_zcta1 <- modzcta_to_zcta %>% mutate(ZCTA = as.character(ZCTA))
modzcta_to_zcta2 <- modzcta_to_zcta1 %>% mutate(MODZCTA = as.character(MODZCTA))
MODZCTAs_in_NYC <- as.character(unique(ZCTA_test_series$MODZCTA))
ZCTAs_in_NYC <- as.character(unique(modzcta_to_zcta$ZCTA))
# download ZIP to Tract crosswalk from HUD
zip_to_tract <- download(
"https://www.huduser.gov/portal/datasets/usps/ZIP_TRACT_062020.xlsx",
"ZIP_TRACT_062020.xlsx",
function(p) suppressWarnings(read_excel(path = p, col_types = c("text", "text", "numeric", "skip", "skip", "skip")))
)
# Block to ZCTA and County crosswalk for NY
ny_xwalk <- download("https://lehd.ces.census.gov/data/lodes/LODES7/ny/ny_xwalk.csv.gz",
"ny_xwalk.csv.gz",
function(p) {
zzf = gzfile(p)
read.csv(zzf) %>%
dplyr::select(cty, tabblk2010, zcta, blklondd, blklatdd) %>%
mutate(tabblk2010 = as.character(tabblk2010),
zcta = as.character(zcta),
cntyfips = as.character(cty)) %>%
dplyr::select(-cty)
}
)
# We have many sources of data, so these just help to combine the various data types
NYC_counties1 <- c("Bronx","Kings","Queens","New York","Richmond")
NYC_counties1_full <- c("Bronx County","Kings County","Queens County","New York County","Richmond County")
NYC_boro_county_match <- tibble(County = c("Bronx","Kings","Queens","New York","Richmond"),
boro = c("Bronx","Brooklyn","Queens","Manhattan","Staten Island"),
full_county = c("Bronx County","Kings County","Queens County","New York County","Richmond County"),
fips = c("36005", "36047", "36081", "36061", "36085"))
# stan model script
BWQS_stan_model <- here("code", "nb_bwqs_cov.stan") Census Data
#### Census Data ####
# function to pull 2010 block population for Queens & Nassau counties
pm(get.qn.blocks <- function(){
nassau_blk_pop <- get_decennial(geography = "block", variables = "P001001",
state = "NY", county = "Nassau", geometry = FALSE)
queens_blk_pop <- get_decennial(geography = "block", variables = "P001001",
state = "NY", county = "Queens", geometry = FALSE)
bind_rows(nassau_blk_pop, queens_blk_pop) %>%
dplyr::select(GEOID, value) %>% rename("pop2010" = "value")
})
# function to pull 2018 ACS data
pm(acs.main <- function(admin_unit = c("zcta", "tract"), state_unit = c(NULL, "NY"), sf_shapes = c(TRUE, FALSE)) {
ACS_Data <- get_acs(geography = admin_unit,
state = state_unit,
geometry = sf_shapes,
variables = c(medincome = "B19013_001",
total_pop1 = "B01003_001",
fpl_100 = "B06012_002",
fpl_100to150 = "B06012_003",
median_rent = "B25031_001",
total_hholds1 = "B22003_001",
hholds_snap = "B22003_002",
over16total_industry1 = "C24050_001",
ag_industry = "C24050_002",
construct_industry = "C24050_003",
manufact_industry = "C24050_004",
wholesaletrade_industry = "C24050_005",
retail_industry = "C24050_006",
transpo_and_utilities_industry = "C24050_007",
information_industry = "C24050_008",
finance_and_realestate_industry = "C24050_009",
science_mngmt_admin_industry = "C24050_010",
edu_health_socasst_industry = "C24050_011",
arts_entertain_rec_accomodate_industry = "C24050_012",
othsvcs_industry = "C24050_013",
publicadmin_industry = "C24050_014",
total_commute1 = "B08301_001",
drove_commute = "B08301_002",
pubtrans_bus_commute = "B08301_011",
pubtrans_subway_commute = "B08301_013",
pubtrans_railroad_commute = "B08301_013",
pubtrans_ferry_commute = "B08301_015",
taxi_commute = "B08301_016",
bicycle_commute = "B08301_018",
walked_commute = "B08301_019",
workhome_commute = "B08301_021",
unemployed = "B23025_005",
under19_noinsurance = "B27010_017",
age19_34_noinsurance = "B27010_033",
age35_64_noinsurance = "B27010_050",
age65plus_noinsurance = "B27010_066",
hisplat_raceethnic = "B03002_012",
nonhispLat_white_raceethnic = "B03002_003",
nonhispLat_black_raceethnic = "B03002_004",
nonhispLat_amerindian_raceethnic = "B03002_005",
nonhispLat_asian_raceethnic = "B03002_006",
age65_plus = "B08101_008"),
year = 2018,
output = "wide",
survey = "acs5")
if(admin_unit=="zcta"){
ACS_Data <- ACS_Data %>% #only pull out the estimates and cleaning variable names
filter(GEOID %in% ZCTAs_in_NYC) %>%
dplyr::select(-NAME) %>%
dplyr::select(GEOID, !ends_with("M")) %>%
rename_at(vars(ends_with("E")), .funs = list(~str_sub(., end = -2)))
}
if(admin_unit=="tract"){
ACS_Data <- ACS_Data %>% #only pull out the estimates and cleaning variable names
filter(substr(GEOID,1,5) %in% NYC_boro_county_match$fips) %>% # Tracts in NYC counties
dplyr::select(-NAME) %>%
dplyr::select(GEOID, !ends_with("M")) %>%
rename_at(vars(ends_with("E")), .funs = list(~str_sub(., end = -2)))
}
return(ACS_Data)
})
# Use 2010 block population to scale Nassau County ZCTA contribution to NYC MODZCTAs
# This is the method used according to data dictionary at NYC Open Data:
# https://data.cityofnewyork.us/Health/Modified-Zip-Code-Tabulation-Areas-MODZCTA-/pri4-ifjk
nassau_zcta_weights <- function(zcta_acs, mz_to_z, blk_to_z){
blk_pop = get.qn.blocks()
modzcta_span = c("11429", "11411", "11004") # MODZCTA with ZCTAs from both Queens & Nassau
border_zcta <- mz_to_z %>% filter(MODZCTA %in% modzcta_span) %>% pull(ZCTA)
# Filter block-zcta crosswalk table to Queens-Nassau ZCTA of interest
blk_to_z <- blk_to_z %>% filter(zcta %in% border_zcta)
# Join population to block-zcta crosswalk
xwalk_pop <- blk_to_z %>% left_join(blk_pop, by = c("tabblk2010" = "GEOID"))
# Summarise 2010 population by ZCTA to calculate proportions inside NYC
zcta_pop_2010 <- xwalk_pop %>%
group_by(zcta) %>%
summarise(z_pop_2010 = sum(pop2010), .groups = "drop_last")
queens_zcta_pop_2010 <- xwalk_pop %>%
filter(cntyfips == "36081") %>%
group_by(zcta) %>%
summarise(queens_z_pop_2010 = sum(pop2010), .groups = "drop_last")
zcta_pop_props <- queens_zcta_pop_2010 %>%
left_join(zcta_pop_2010, by = "zcta") %>%
mutate(in_NYC_prop = queens_z_pop_2010/z_pop_2010)
zcta_weights <- zcta_acs %>% dplyr::select(GEOID) %>%
left_join(zcta_pop_props, by = c("GEOID" = "zcta")) %>%
dplyr::select(GEOID, in_NYC_prop) %>%
mutate(in_NYC_prop =
case_when(is.na(in_NYC_prop) ~ 1,
TRUE ~ in_NYC_prop))
# Apply weights for all Census variables except median vars
zcta_acs <- zcta_acs %>% left_join(zcta_weights, by = "GEOID")
varvec <- 1:ncol(zcta_acs)
varvec <- varvec[-grep("GEOID|in_NYC_prop|medincome|median_rent", names(zcta_acs))]
zcta_acs <- zcta_acs %>% mutate_at(vars(all_of(varvec)), ~ . * in_NYC_prop)
}
# function to clean ACS data
clean_acs_data_and_derive_vars <- function(df, admin_unit = c("zcta", "tract")){
if(admin_unit=="zcta"){
ACS_Data1a <- df %>%
left_join(., modzcta_to_zcta1, by = c("GEOID" = "ZCTA")) %>%
group_by(MODZCTA) %>%
summarise_at(vars(medincome, median_rent), ~weighted.mean(., total_pop1, na.rm = T)) %>%
rename(zcta = "MODZCTA")
ACS_Data1b <- df %>%
left_join(., modzcta_to_zcta1, by = c("GEOID" = "ZCTA")) %>%
group_by(MODZCTA) %>%
summarise_at(vars(total_pop1:fpl_100to150, total_hholds1:age65_plus), ~sum(.)) %>%
mutate_at(vars(ends_with("_commute")), ~round((./total_commute1)*100, 2)) %>% #proportion of people relying on a given mode of transit
mutate_at(vars(ends_with("_raceethnic")), ~round((./total_pop1)*100, 2)) %>% #proportion of ppl reporting a given race/ethncity
mutate(not_insured = round(((under19_noinsurance + age19_34_noinsurance + age35_64_noinsurance + age65plus_noinsurance) / total_pop1)*100, 2), #proportion uninsured
#snap_hholds = round((hholds_snap/total_hholds1)*100, 2), #proportion relying on SNAP
#fpl_150 = round(((fpl_100+fpl_100to150)/total_pop1)*100, 2), #proportion 150% or less of FPL
unemployed = round((unemployed/over16total_industry1)*100, 2), #proportion unemployed
not_quarantined_jobs = round(((ag_industry+(construct_industry*.25)+wholesaletrade_industry+ #an estimate of who is still leaving the house for work based on industry
(edu_health_socasst_industry*.5)+transpo_and_utilities_industry)/over16total_industry1)*100, 2)) %>%
dplyr::select(-ends_with("_noinsurance"), -fpl_100, -fpl_100to150, -ends_with("_industry"), -hholds_snap) %>%
rename(zcta = "MODZCTA")
ACS_Data2 <- left_join(ACS_Data1a, ACS_Data1b, by = "zcta") %>%
mutate(zcta = as.character(zcta))
} else{
ACS_Data2 <- df %>%
mutate_at(vars(ends_with("_commute")), ~round((./total_commute1)*100, 2)) %>% #proportion of people relying on a givenmode of transit
mutate_at(vars(ends_with("_raceethnic")), ~round((./total_pop1)*100, 2)) %>% #proportion of ppl reporting a given race/ethncity
mutate(not_insured = round(((under19_noinsurance + age19_34_noinsurance + age35_64_noinsurance + age65plus_noinsurance) / total_pop1)*100, 2), #proportion uninsured
unemployed = round((unemployed/over16total_industry1)*100, 2), #proportion unemployed
not_quarantined_jobs = round(((ag_industry+(construct_industry*.25)+wholesaletrade_industry+ #an estimate of who is still leaving the house for work based on industry
(edu_health_socasst_industry*.5)+transpo_and_utilities_industry)/over16total_industry1)*100, 2)) %>%
dplyr::select(-ends_with("_noinsurance"), -ends_with("_industry"),-fpl_100, -fpl_100to150,-hholds_snap)
}
return(ACS_Data2)
}
# Functions to pull mode of transportation for our approximate of essential workers
pm(fst = T,
get_essential_acs <- function(admin_unit, state_unit) {
get_acs(geography = admin_unit, #pull down the relevant categories
state = state_unit,
variables = c(ag_car1_commute = "B08126_017",
ag_pubtrans_commute = "B08126_047",
construct_car1_commute ="B08126_018",
construct_pubtrans_commute = "B08126_048",
wholesale_car1_commute = "B08126_020",
wholesale_pubtrans_commute = "B08126_050",
transpo_car1_commute = "B08126_022",
transpo_pubtrans_commute = "B08126_052",
ed_hlthcare_car1_commute = "B08126_026",
ed_hlthcare_pubtrans_commute = "B08126_056"),
year = 2018,
output = "wide",
survey = "acs5")
})
acs.essential <- function(admin_unit, zcta_pop = NULL, state_unit = NULL) {
if(!admin_unit %in% c("zcta", "tract")) stop("admin_unit must be either 'zcta' or 'tract'")
if(admin_unit == "tract" & is.null(state_unit)) stop("state_unit must be set to download tracts")
if(!is.null(state_unit)) if(state_unit != "NY") stop("state_unit must be either NULL or 'NY'")
ACS_EssentialWrkr_Commute <- get_essential_acs(admin_unit = admin_unit, state_unit = state_unit)
if(admin_unit == "zcta"){
if(is.null(zcta_pop)) stop("zcta_pop must be set for scaling MODZCTA on Queens/Nassau boundary")
ACS_Essential_worker_estimates <- ACS_EssentialWrkr_Commute %>% #clean data and aggregate
dplyr::select(-ends_with("M"), -NAME) %>%
filter(GEOID %in% ZCTAs_in_NYC) %>%
# scale ZCTAs by proportion of population in NYC
right_join(zcta_pop %>% dplyr::select("GEOID", "in_NYC_prop"), by = "GEOID") %>%
mutate_at(vars(2:11), ~ . * in_NYC_prop) %>%
# summarize ZCTA to MODZCTA
left_join(., modzcta_to_zcta1, by = c("GEOID" = "ZCTA")) %>%
group_by(MODZCTA) %>%
summarise_at(vars(2:11),
~ sum(., na.rm = T)) %>%
rename(zcta = "MODZCTA") %>%
mutate_at(vars(starts_with("ed_hlthcare")), ~ round(. / 2), 0) %>% #maintain same proportions as estimated nonquarintined jobs
mutate_at(vars(starts_with("construct")), ~ round(. / 4), 0) %>%
mutate(
essentialworker_drove = rowSums(dplyr::select(., contains("car1_commute"))),
essentialworker_pubtrans = rowSums(dplyr::select(., contains("pubtrans")))) %>%
dplyr::select(zcta, essentialworker_drove, essentialworker_pubtrans) %>%
mutate(zcta = as.character(zcta))
}
else { # tracts
ACS_Essential_worker_estimates <- ACS_EssentialWrkr_Commute %>% #clean data and aggregate
dplyr::select(-ends_with("M"), -NAME) %>%
filter(substr(GEOID,1,5) %in% NYC_boro_county_match$fips) %>% # Tracts in NYC counties
mutate_at(vars(starts_with("ed_hlthcare")), ~round(./2), 0) %>% #maintain same proportions as estimated nonquarintined jobs
mutate_at(vars(starts_with("construct")), ~round(./4), 0) %>%
mutate(essentialworker_drove = rowSums(dplyr::select(., contains("car1_commute"))),
essentialworker_pubtrans = rowSums(dplyr::select(., contains("pubtrans")))) %>%
dplyr::select(GEOID, essentialworker_drove, essentialworker_pubtrans)
}
return(ACS_Essential_worker_estimates)
}
# ZCTA CENSUS DATA
options(tigris_use_cache = TRUE)
ACS_Data1 <- as.data.frame(acs.main("zcta", NULL, FALSE)) #download the zcta data## Getting data from the 2014-2018 5-year ACSACS_Data_scaled <- nassau_zcta_weights(ACS_Data1, modzcta_to_zcta2, ny_xwalk)## Getting data from the 2010 decennial Census## Getting data from the 2010 decennial CensusACS_Data2 <- clean_acs_data_and_derive_vars(ACS_Data_scaled, "zcta")
ACS_EssentialWrkr_Commute1 = as.data.frame(acs.essential("zcta", zcta_pop = ACS_Data_scaled, state_unit = NULL))## Getting data from the 2014-2018 5-year ACSprint(paste("The 2018 5-year ACS population range in NYC MODZCTAs is:", paste(range(ACS_Data2$total_pop1), collapse = "-")))## [1] "The 2018 5-year ACS population range in NYC MODZCTAs is: 3028-112425"# TRACT CENSUS DATA
acs_tracts <- acs.main("tract", "NY", TRUE)## Getting data from the 2014-2018 5-year ACSacs_tracts2 <- clean_acs_data_and_derive_vars(acs_tracts, "tract")
acs_tracts_commute1 = as.data.frame(acs.essential("tract", state_unit = "NY"))## Getting data from the 2014-2018 5-year ACSt_census = Sys.time()Grocery stores per area
#### Grocery stores per area ####
non_supermarket_strings <- c("DELI|TOBACCO|GAS|CANDY|7 ELEVEN|7-ELEVEN|LIQUOR|ALCOHOL|BAKERY|CHOCOLATE|DUANE READE|WALGREENS|CVS|RITE AID|RAVIOLI|WINERY|WINE|BEER|CAFE|COFFEE")
food_retail_filtered <- food_retail %>%
filter(County %in% NYC_boro_county_match$County) %>%
filter(str_detect(`Establishment Type`, "J") & str_detect(`Establishment Type`, "A") & str_detect(`Establishment Type`, "C") &
!str_detect(`Establishment Type`, "H")) %>%
filter(!str_detect(`Entity Name`, non_supermarket_strings) & !str_detect(`DBA Name`, non_supermarket_strings)) %>%
filter(`Square Footage`>=4500) %>%
mutate(zcta = as.character(str_extract(Location, "[:digit:]{5}"))) %>%
mutate(Address = case_when(
# some locations geocode better when address includes city name
`License Number` %in% c("638599", "712410", "706967", "710078") ~
paste(paste(`Street Number`, `Street Name`), City, State, zcta, sep = ","),
# but most geocode better without it, see limitation #4 at: http://gis.ny.gov/gisdata/inventories/details.cfm?DSID=1278
TRUE ~
paste(paste(`Street Number`, `Street Name`), State, zcta, sep = ",") )
)
# Geocode grocers, using a cached version if available to make analysis reproducible
# The geocoding service may be updated in the future and give different results
cached_grocers = file.path(data.root, "grocers_geocode_2020-11-09.csv")
if(file.exists(cached_grocers) & use_repo_data){
gctable <- read.csv(cached_grocers)
failed = which(gctable$score == 0)
message("Loaded cached geocoded grocers: ", nrow(gctable)-length(failed), "/", nrow(gctable), " have coordinates.") # 997/1037
if(nrow(gctable) != nrow(food_retail_filtered)) warning("Cached geocoded table has different row count than non-geocoded table")
} else {
# locations are returned in crs=26918, UTM 18N NAD83
api = "https://gisservices.its.ny.gov/arcgis/rest/services/Locators/Street_and_Address_Composite/GeocodeServer/findAddressCandidates?f=json&maxLocations=1&SingleLine="
message("Geocoding ", nrow(food_retail_filtered), " grocers with NY State ITS geocoder...")
t1 = Sys.time()
res = lapply(food_retail_filtered$Address, function(addr) {
GET(url = paste0(api, URLencode(addr))) })
tdiff = Sys.time() - t1
# extract results
geocodes = lapply(res, function(page) fromJSON(rawToChar(page$content), flatten = TRUE)$candidates)
failed = which(sapply(geocodes,class) != "data.frame")
geocodes[failed] <- lapply(1:length(failed), function(void){
data.frame(address = NA_character_, score = 0, location.x = NA_real_, location.y = NA_real_)})
gctable = bind_rows(geocodes)
message("Geocoded ", nrow(food_retail_filtered)-length(failed), "/", nrow(food_retail_filtered),
" grocers in ", round(as.numeric(tdiff),1), " ", attributes(tdiff)$units)
write.csv(gctable, paste0(data.root, "/grocers_geocode_", Sys.Date(), ".csv"))
}## Loaded cached geocoded grocers: 997/1037 have coordinates.# Count grocers by tract
gctable = filter(gctable, score > 0)
grocerSF = st_as_sf(gctable, coords = c("location.x", "location.y"), crs = 26918) %>% st_transform(crs = 2263)
tractSF = acs_tracts2[, "GEOID"] %>% st_transform(., crs = 2263)
tract_grocers = suppressWarnings(st_intersection(tractSF, grocerSF)) %>%
st_set_geometry(., NULL) %>%
group_by(GEOID) %>%
summarise(grocers = n_distinct(`address`))
#nrow(tract_grocers) # 749
#sum(tract_grocers$grocers) # 993
# Count grocers by ZCTA
zctaSF <- MODZCTA_NYC_shp %>% dplyr::select(modzcta, geometry) %>% st_transform(crs = 2263)
zcta_grocers <- suppressWarnings(st_intersection(zctaSF, grocerSF)) %>%
st_set_geometry(., NULL) %>%
group_by(modzcta) %>%
summarise(grocers = n_distinct(`address`))
# sum(zcta_grocers$grocers) # 993
# nrow(zcta_grocers) # 172
# range(zcta_grocers$grocers) # 1, 21Subway station locations
#### Subway station locations ####
SubwayStation_shp <- as_tibble(turnstile()$stations) %>%
st_as_sf(., coords = c("lon", "lat"), crs = 4269) %>%
st_transform(., crs = 2263) %>%
filter(!str_detect(ca, "PTH")) #removing New Jersey PATH stationsResidential area
#### Residential area ####
Pluto_ResOnly <- Pluto %>%
filter(landuse>="01" & landuse<="04") %>%
mutate(base_bbl = as.character(bbl)) %>%
dplyr::select(-bbl)
ResBBLs <- as.character(Pluto_ResOnly$base_bbl)
# zcta-level residential area
Res_Bldg_Footprints <- Bldg_Footprints %>%
st_set_geometry(., NULL) %>%
mutate(base_bbl = as.character(base_bbl)) %>%
filter(base_bbl %in% ResBBLs &
feat_code == "2100") %>%
mutate(bldg_volume = shape_area * heightroof) %>%
left_join(., Pluto_ResOnly, by = "base_bbl") %>%
mutate(bldg_volume = if_else(is.na(bldg_volume), shape_area*numfloors*10, bldg_volume),
res_volume = (bldg_volume/unitstotal)*unitsres,
zcta = as.character(zipcode)) %>%
group_by(zcta) %>%
summarise(total_res_volume_zcta = sum(res_volume, na.rm = TRUE))
# tract-level residential area
Res_Bldg_Footprints2 <- Bldg_Footprints %>%
st_transform(crs = 2263) %>%
suppressWarnings(st_centroid(of_largest_polygon = TRUE)) %>%
mutate(base_bbl = as.character(base_bbl)) %>%
filter(base_bbl %in% ResBBLs &
feat_code == "2100") %>%
mutate(bldg_volume = shape_area * heightroof) %>%
left_join(., Pluto_ResOnly, by = "base_bbl") %>%
mutate(bldg_volume = if_else(is.na(bldg_volume), shape_area*numfloors*10, bldg_volume),
res_volume = (bldg_volume/unitstotal)*unitsres)
pm(get.tract.res <- function(res, tracts) st_intersection(res, tracts)) # takes a few minutes
res_bldg_tract <- get.tract.res(Res_Bldg_Footprints2, tractSF)## Warning: attribute variables are assumed to be spatially constant throughout all
## geometriesres_bldg_tract_sum <- st_set_geometry(res_bldg_tract, NULL) %>%
group_by(GEOID) %>%
summarise(total_res_volume_tract = sum(res_volume, na.rm = TRUE))
#nrow(res_bldg_tract_sum) # 2132COVID Tests
#### COVID Tests ####
MODZCTA_NYC_shp1 <- MODZCTA_NYC_shp %>%
dplyr::select(modzcta, geometry) %>%
rename("zcta" = "modzcta")
May7_tests <- ZCTA_test_series %>%
filter(date=="2020-05-07") %>%
mutate(zcta = as.character(MODZCTA)) %>%
rename(total_tests = "Total") %>%
dplyr::select(zcta, date, Positive, total_tests)
ZCTA_by_boro1 <- ZCTA_by_boro %>%
mutate(boro = as.character(boro),
zcta = as.character(zip)) %>%
dplyr::select(-zip) %>%
bind_rows(.,
tibble(boro = as.character(c("Manhattan", "Manhattan" ,"Queens")), #correcting nas
zcta = as.character(c("10069", "10282", "11109"))))
# get water mask for maps
source(here("code/water_mask.R"))## although coordinates are longitude/latitude, st_intersects assumes that they are planar# Figure 3B - Map of Tests by MODZCTA
theme_smallmaps <- theme(legend.title = element_text(face = "bold", size = 9),
plot.title = element_text(size = 9.5),
panel.background = element_rect(fill = "#cccccc"),
panel.grid = element_blank(),
legend.background = element_rect(fill = "transparent"),
legend.position = c(0.22, 0.80),
legend.text = element_text(size = 8.5),
plot.margin = unit(c(4,0,4,0), "pt"),
legend.key.size = unit(1.1, "lines"),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.ticks.length = unit(0, "pt"),
axis.title = element_blank())
fig3b <- MODZCTA_NYC_shp1 %>%
left_join(., May7_tests, by = "zcta") %>%
left_join(., ACS_Data2, by = "zcta") %>%
filter(zcta != "99999") %>%
mutate(pos_per_100000 = (Positive/total_pop1)*100000) %>%
ggplot() +
geom_sf(data = basemap_water, fill = "white", lwd = 0) +
geom_sf(aes(fill = pos_per_100000), lwd = 0.2)+
scalebar(MODZCTA_NYC_shp1, dist = 5, dist_unit = "km",
transform = TRUE, model = "WGS84",
st.size = 2.8, height = 0.015, border.size = 0.5,
anchor = c(x = -73.71, y = 40.51)) +
labs(fill = "Positives per 100,000") +
ggtitle("Cumulative positive COVID tests by zip code (May 7, 2020)") +
scale_fill_gradientn(colours=brewer_pal("YlGnBu", type = "seq")(7)) +
coord_sf(crs = st_crs(MODZCTA_NYC_shp1),
xlim = c(plot_bounds$xmin, plot_bounds$xmax),
ylim = c(plot_bounds$ymin, plot_bounds$ymax),
expand = FALSE) +
theme_bw(base_size = 6) +
theme_smallmapsCreate data frames of all above information
#### Create data frames of all above information ####
ZCTA_ACS_COVID_shp <- MODZCTA_NYC_shp1 %>%
st_transform(., crs = 2263) %>%
dplyr::select(zcta, geometry) %>%
left_join(., ACS_Data2, by = "zcta") %>%
left_join(., May7_tests, by = "zcta") %>%
left_join(., Res_Bldg_Footprints, by = "zcta") %>%
left_join(., ACS_EssentialWrkr_Commute1, by = "zcta") %>%
left_join(., zcta_grocers, by = c("zcta" = "modzcta")) %>%
mutate(pop_density = as.numeric(total_pop1/st_area(geometry)),
avg_hhold_size = round((total_pop1/total_hholds1), 2),
pos_per_100000 = (Positive/total_pop1)*100000,
testing_ratio = (total_tests/total_pop1),
res_vol_zctadensity = as.numeric(total_res_volume_zcta/st_area(geometry)),
res_vol_popdensity = as.numeric(total_pop1/total_res_volume_zcta),
pubtrans_ferrysubway_commute = pubtrans_subway_commute + pubtrans_ferry_commute,
grocers = replace_na(grocers, 0),
grocers_per_1000 = (grocers/total_pop1)*1000,
pos_per_100000 = round(pos_per_100000, 0),
valid_var = "0",
one_over_grocers_per_1000 = if_else(is.infinite(1/grocers_per_1000), (1/.0293)+1, 1/grocers_per_1000), #max + 1 for zero values
didnot_workhome_commute = 100 - workhome_commute,
one_over_medincome = 1/medincome) %>%
dplyr::select(-pubtrans_subway_commute, -pubtrans_ferry_commute) %>%
left_join(., ZCTA_by_boro1, by = "zcta") %>%
mutate_at(vars(starts_with("essentialworker_")), ~round((./over16total_industry1)*100, 2)) %>%
filter(zcta != "99999") #remove na
ZCTA_ACS_COVID <- ZCTA_ACS_COVID_shp %>%
st_set_geometry(., NULL) #remove geometry
tract_vars <- tractSF %>% # uses local CRS
left_join(., st_set_geometry(acs_tracts2, NULL), by = "GEOID") %>%
left_join(., acs_tracts_commute1, by = "GEOID") %>%
left_join(., res_bldg_tract_sum, by = "GEOID") %>%
left_join(., tract_grocers, by = "GEOID") %>%
mutate(pop_density = as.numeric(total_pop1/st_area(geometry)),
avg_hhold_size = round((total_pop1/total_hholds1), 2),
res_vol_tractdensity = as.numeric(total_res_volume_tract/st_area(geometry)),
res_vol_popdensity = as.numeric(total_pop1/total_res_volume_tract),
pubtrans_ferrysubway_commute = pubtrans_subway_commute + pubtrans_ferry_commute,
grocers = replace_na(grocers, 0),
grocers_per_1000 = (grocers/total_pop1)*1000,
valid_var = "0",
didnot_workhome_commute = 100 - workhome_commute,
one_over_grocers_per_1000 = if_else(is.infinite(1/grocers_per_1000), (1/.0293)+1, 1/grocers_per_1000),
one_over_medincome = 1/medincome) %>%
dplyr::select(-pubtrans_subway_commute, -pubtrans_ferry_commute) %>%
mutate_at(vars(starts_with("essentialworker_")), ~round((./over16total_industry1)*100, 2))Tract to ZCTA weighted assignment
#### Tract to ZCTA weighted assignment ####
# Shares of residential addresses in ZCTAs by Tract are later used to select
# representative tract-level SES values at a specified location on the ECDF
modzcta_to_zcta_chr <- data.frame(ZCTA = as.character(modzcta_to_zcta$ZCTA),
MODZCTA = as.character(modzcta_to_zcta$MODZCTA))
# Calculate the proportion of population each combined ZCTA contributes to MODZCTA
modzcta_to_zcta_pop <- modzcta_to_zcta_chr %>%
left_join(ACS_Data_scaled[, c("GEOID", "total_pop1")], by = c("ZCTA" = "GEOID")) %>%
group_by(MODZCTA) %>%
mutate(sum_pop = sum(total_pop1)) %>%
ungroup() %>%
mutate(pop_prop = total_pop1/sum_pop) %>%
arrange(MODZCTA)
# Sum ZCTA->Tract RES_RATIO weights when multiple ZCTA combine to a single MODZCTA
modzcta_to_tract <- modzcta_to_zcta_pop %>%
filter(MODZCTA != "99999") %>%
# Weights of each ZCTA are scaled by their proportion of MODZCTA population
left_join(zip_to_tract, by = c("ZCTA" = "ZIP")) %>%
mutate(scaled_res_ratio = RES_RATIO * pop_prop) %>%
# note: MODZCTA 10001 & 10007 contain a ZCTA with a small share of population
# but a zero share of RES_RATIO. This makes their sum of scaled_res_ratio by
# MODZCTA less than 1, but this will not effect the weighted median or 3Q values
# by tract.
group_by(MODZCTA, TRACT) %>%
dplyr::summarize(SUM_RES_RATIO = sum(scaled_res_ratio), .groups = "drop") %>%
# Weights are then multiplied by the population density of each tract
left_join(tract_vars, by = c("TRACT" = "GEOID")) %>%
dplyr::select(MODZCTA, TRACT, SUM_RES_RATIO, total_pop1, total_hholds1) %>%
filter(total_hholds1 > 0) %>%
mutate(tract_popdens = total_pop1/total_hholds1) %>%
mutate(W_RES_RATIO = SUM_RES_RATIO * tract_popdens) %>%
dplyr::select(MODZCTA, TRACT, W_RES_RATIO)
modzcta_to_tract## # A tibble: 3,020 x 3
## MODZCTA TRACT W_RES_RATIO
## <chr> <chr> <dbl>
## 1 10001 36061005600 0.0307
## 2 10001 36061005800 0.0648
## 3 10001 36061007100 0
## 4 10001 36061007600 0.250
## 5 10001 36061008400 0.00422
## 6 10001 36061009100 0.154
## 7 10001 36061009300 0.145
## 8 10001 36061009500 0.231
## 9 10001 36061009700 0.250
## 10 10001 36061009900 0.253
## # … with 3,010 more rows# length(unique(modzcta_to_tract$MODZCTA)) # 177
# length(unique(modzcta_to_tract$TRACT)) # 2111
# length(unique(tractSF$GEOID)) # 2167
# check res_ratio against tracts with no population
tract_modzcta_pop <- modzcta_to_tract %>%
left_join(acs_tracts2, by = c("TRACT" = "GEOID")) %>%
dplyr::select(MODZCTA, TRACT, total_pop1, W_RES_RATIO)
# checking for tracts with no population but res_ratio > 0
# (not possible now that the weighted ratio uses population)
#tract_modzcta_pop %>% filter(total_pop1 == 0 & W_RES_RATIO > 0) %>% nrow() # 0
# check to make sure no MODZCTA have all zeroes for all tract res_ratios
all_zero_ratio <- tract_modzcta_pop %>% group_by(MODZCTA) %>% filter(all(W_RES_RATIO == 0))
if(nrow(all_zero_ratio) > 0){
warning("The following MODZCTA have all zeroes for weighted tract residents:\n ",
paste(all_zero_ratio$MODZCTA, collapse = ","))
}
modzcta_to_tract2 <- dplyr::select(tract_modzcta_pop, -total_pop1)
### Preparation done --- Now for the Analysis ###Part 1: Creation of BWQS Neighborhood Infection Risk Scores
#### Part 1: Creation of BWQS Neighborhood Infection Risk Scores ####
# Step 1: Create univariate scatterplots with negative binomial linear smoothers
# to make sure direction of associations are consistent for all variables
# Visualize the non-linear relationship between infection rate and test_ratio (not adjusted for social inequity)
# infections versus testing ratio with smoothers; this covariate is not quantiled so the shape of association is important
ggplot(ZCTA_ACS_COVID, aes(x = testing_ratio, y = pos_per_100000)) + geom_point() +
ylab("infections per 100,000") +
geom_smooth(color = "red", formula = y ~ x, method = "glm.nb") +
stat_smooth(color = "green", method = "gam", formula = y ~ s(x), method.args = list(family = "nb"), se = F) +
geom_smooth(method = "glm.nb", formula = y ~ splines::ns(x, 3), se = FALSE) + theme_bw()
# 10 social variables
plotsocial <- ggplot(ZCTA_ACS_COVID, aes(y = pos_per_100000)) +
geom_point() + geom_smooth(color = "red", formula = y ~ x, method = "glm.nb") + ylab("infections per 100,000") + theme_bw()
plotsocial %+% aes(x = one_over_medincome)
plotsocial %+% aes(x = not_insured)
plotsocial %+% aes(x = one_over_grocers_per_1000)
plotsocial %+% aes(x = unemployed)
plotsocial %+% aes(x = res_vol_popdensity)
plotsocial %+% aes(x = didnot_workhome_commute)
plotsocial %+% aes(x = not_quarantined_jobs)
plotsocial %+% aes(x = avg_hhold_size)
plotsocial %+% aes(x = essentialworker_drove)
plotsocial %+% aes(x = essentialworker_pubtrans)
SES_vars <- names(ZCTA_ACS_COVID %>% dplyr::select(one_over_medincome, not_insured, one_over_grocers_per_1000, unemployed,
not_quarantined_jobs, essentialworker_pubtrans, essentialworker_drove,
didnot_workhome_commute, res_vol_popdensity, avg_hhold_size))
# Step 2a: Examine relationships between explanatory variables to make sure nothing >0.9 correlation, as this could bias BWQS
Cors_SESVars <- cor(x = ZCTA_ACS_COVID %>% dplyr::select(all_of(SES_vars)), method = "kendall")
Cors_SESVars1 <- as.data.frame(Cors_SESVars)
Cors_SESVars1$var1 <- row.names(Cors_SESVars1)
Cors_SESVars2 <- gather(data = Cors_SESVars1, key = "var2", value = "correlation", -var1) %>%
filter(var1 != var2)
Cors_SESVars2 %>% arrange(desc(correlation)) %>% distinct(correlation, .keep_all = T) ## var1 var2 correlation
## 1 essentialworker_drove not_quarantined_jobs 0.6121033
## 2 essentialworker_pubtrans one_over_medincome 0.5432472
## 3 unemployed one_over_medincome 0.5341512
## 4 not_insured one_over_medincome 0.5270258
## 5 essentialworker_pubtrans unemployed 0.5243542
## 6 didnot_workhome_commute essentialworker_drove 0.5086620
## 7 avg_hhold_size not_quarantined_jobs 0.4787113
## 8 didnot_workhome_commute not_quarantined_jobs 0.4763680
## 9 avg_hhold_size didnot_workhome_commute 0.4598744
## 10 res_vol_popdensity one_over_medincome 0.4549307
## 11 res_vol_popdensity essentialworker_pubtrans 0.4507468
## 12 res_vol_popdensity not_insured 0.4391560
## 13 avg_hhold_size essentialworker_drove 0.4198843
## [ reached 'max' / getOption("max.print") -- omitted 32 rows ]if(export.figs) {
png(filename = file.path(fig.path, paste0("sfig1_", Sys.Date(), ".png")), width = 96*7, height = 96*7)
corrplot(Cors_SESVars, p.mat = Cors_SESVars, insig = "p-value", type = "lower",
sig.level = -1, tl.col = "black", tl.srt = 45)
dev.off()
}## png
## 2
# Step 2b: Examine Univariable kendall associations for all selected variables with the outcome
bind_cols(Variables = SES_vars,
ZCTA_ACS_COVID %>%
dplyr::select(all_of(SES_vars), pos_per_100000) %>%
summarise_at(vars(all_of(SES_vars)), list(~cor(., pos_per_100000, method = "kendall"))) %>%
t() %>%
as_tibble(),
ZCTA_ACS_COVID %>%
dplyr::select(all_of(SES_vars), pos_per_100000) %>%
summarise_at(vars(all_of(SES_vars)),
list(~cor.test(., pos_per_100000, method = "kendall")$p.value)) %>%
t() %>%
as_tibble()) %>%
mutate(`Correlation (Tau)` = round(V1...2, 3),
`p value` = as.character(ifelse(V1...3 < 0.0001, "< 0.0001", round(V1...3, 3))),) %>%
dplyr::select(-V1...2, -V1...3) ## Warning: The `x` argument of `as_tibble.matrix()` must have column names if `.name_repair` is omitted as of tibble 2.0.0.
## Using compatibility `.name_repair`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.## New names:
## * V1 -> V1...2
## * V1 -> V1...3## # A tibble: 10 x 3
## Variables `Correlation (Tau)` `p value`
## <chr> <dbl> <chr>
## 1 one_over_medincome 0.322 < 0.0001
## 2 not_insured 0.236 < 0.0001
## 3 one_over_grocers_per_1000 0.268 < 0.0001
## 4 unemployed 0.294 < 0.0001
## 5 not_quarantined_jobs 0.542 < 0.0001
## 6 essentialworker_pubtrans 0.175 0.001
## 7 essentialworker_drove 0.483 < 0.0001
## 8 didnot_workhome_commute 0.453 < 0.0001
## 9 res_vol_popdensity 0.106 0.035
## 10 avg_hhold_size 0.463 < 0.0001# Step 3: Prepare data for BWQS and pass to stan for model fitting
pm(fit_BWQS_model <- function(data_list, stan_model_path){
fitted_model <- stan(
file = stan_model_path,
data = data_list,
chains = 1,
warmup = 2500,
iter = 20000,
thin = 10,
refresh = 0,
algorithm = "NUTS",
seed = 1234,
control = list(max_treedepth = 20,
adapt_delta = 0.999999999999999))
return(fitted_model)
})
Compute_Bayes_R2 <- function(fit) {
y_pred = exp(extract(fit,"eta")$eta)
var_fit = apply(y_pred, 1, var)
mean_fit = apply(y_pred, 1, mean)
phi = extract(fit,"phi")$phi
var_res = mean_fit + (mean_fit^2)/phi
r2 = var_fit / (var_fit + var_res)
return(list(meanr2 = mean(r2),
medianr2 = median(r2),
lowerr2 = quantile(r2, .025),
upperr2 = quantile(r2, .975)))
}
prep_BWQS_data <- function(df, ses_varnames){
y = as.numeric(df$pos_per_100000)
X <- df %>%
dplyr::select(all_of(ses_varnames))
K <- ns(df$testing_ratio, df = 3)
for (vname in ses_varnames)
X[[vname]] <- ecdf(X[[vname]])(X[[vname]]) * 10
data <-as.data.frame(cbind(y,X)) # Aggregate data in a data.frame
data_list = list(N = NROW(data),
N_new = NROW(data),
C = NCOL(X),
K = NCOL(K),
XC = cbind(as.matrix(X)),
XK = cbind(K),
XC_new = cbind(as.matrix(X)),
XK_new = cbind(K),
Dalp = rep(1,length(ses_varnames)),
y = as.vector(data$y))
return(list(data_list = data_list, X = X, K = K))
}
m1data <- prep_BWQS_data(ZCTA_ACS_COVID, SES_vars)
t_dataprep = Sys.time()
# fit our primary model -- negative binomial
m1 <- fit_BWQS_model(m1data$data_list, BWQS_stan_model)
t_m1 = Sys.time()
# print m1 for model diagnostics (n_eff as % of 1750, Rhat, trace plots, acf)
# m1
min(summary(m1)$summary[,"n_eff"]/1750) # effective sample size is at worst 78%## [1] 0.6186883traceplot(m1, pars = c("beta1", "W"))
stan_ac(m1, pars = c("beta1", "W"))## Warning: Ignoring unknown parameters: fun.y## No summary function supplied, defaulting to `mean_se()`## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
## No summary function supplied, defaulting to `mean_se()`
extract_waic(m1)## $waic
## waic
## 2369.804
##
## $elpd_waic
## elpd_waic
## -1184.902
##
## $p_waic
## p_waic
## 14.69814
##
## $elpd_loo
## elpd_loo
## -1184.963
##
## $p_loo
## p_loo
## 14.75942round(Compute_Bayes_R2(m1)$meanr2, 2)## [1] 0.93round(Compute_Bayes_R2(m1)$lowerr2, 2)## 2.5%
## 0.92round(Compute_Bayes_R2(m1)$upperr2, 2)## 97.5%
## 0.95# residual analysis with DHARMa
residuals_qq <- createDHARMa(simulatedResponse = t(extract(m1,"y_new")$y_new), observedResponse = m1data$data_list$y)## No fitted predicted response provided, using the mean of the simulations#plotQQunif(residuals_qq, testOutliers = F, testDispersion = F) #base graphicsks_unif <- testUniformity(residuals_qq)residuals_qq_unif <- gap::qqunif(residuals_qq$scaledResiduals,pch=2,bty="n", logscale = F, col = "black", cex = 0.6, main = "QQ plot residuals", cex.main = 1)sfig3_qq <- ggplot(data.frame(x = residuals_qq_unif$x, y = residuals_qq_unif$y)) +
geom_abline(linetype = "dashed", color = "blue") +
scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
scale_y_continuous(limits = c(0, 1), expand = c(0, 0)) +
geom_point(aes(x, y), size = 0.8) +
annotate("text", x = 0.02, y = 0.9, hjust = 0, size = 3.2,
label = paste0("Kolmogorov-Smirnov test for\ncomparison of distributions: p=",
round(ks_unif$p.value, 2))) +
coord_equal(ratio = 1) +
labs(x = "Expected", y = "Observed") +
theme_bw() + theme(plot.margin = unit(c(5.5, 11, 5.5, 5.5), "points"))
if(export.figs) ggsave(filename = file.path(fig.path, paste0("sfig3_", Sys.Date(), ".png")), width = 3.5, height = 3.4)
# examine parameter estimates
exp(mean(extract(m1, "beta1")$beta1))## [1] 1.076975vars = c("phi", "beta0", "beta1", paste0("delta", 1:3), SES_vars)
parameters_to_drop <- c("phi", paste0("delta", 1:3), "beta0", "beta1")
number_of_coefficients <- length(vars)
BWQS_params <- bind_cols(as_tibble(summary(m1)$summary[1:number_of_coefficients,c(1,4,6,8:10)]), label = vars)
BWQS_params %>% filter(label == "beta1") %>% mutate_at(vars(1:3), ~exp(.))## # A tibble: 1 x 7
## mean `2.5%` `50%` `97.5%` n_eff Rhat label
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 1.08 1.06 1.08 0.0855 1631. 1.00 beta1BWQS_fits <- BWQS_params %>%
rename(lower = "2.5%",
upper = "97.5%") %>%
filter(!label %in% parameters_to_drop) %>%
arrange(desc(mean)) %>%
mutate(group = factor(if_else(label == "one_over_medincome"|label =="not_insured"|label == "unemployed", "Finances &\nAccess to care",
if_else(label == "one_over_grocers_per_1000", "Food\nAccess",
if_else(str_detect(label, "essential")|label == "not_quarantined_jobs"|label=="didnot_workhome_commute", "Commuting and\nEssential Work",
if_else(label == "avg_hhold_size"|label == "res_vol_popdensity", "Population\nDensity", "Unmatched")))),
levels = c("Commuting and\nEssential Work", "Finances &\nAccess to care", "Population\nDensity", "Food\nAccess")))
labels1 <- c("one_over_medincome" = "1/\nMedian income",
"not_insured" = "Uninsured",
"unemployed" = "Unemployed",
"one_over_grocers_per_1000" = "1/\nGrocers per 1000",
"not_quarantined_jobs" = "Essential Workers",
"essentialworker_pubtrans" = "Essential Worker:\nPublic Transit",
"essentialworker_drove" = "Essential Worker:\nDriving Commute",
"didnot_workhome_commute" = "1/\nWork from home",
"res_vol_popdensity" = "Population Density\nby Residential Volume",
"avg_hhold_size" = "Average people\nper household")
labels2 <- c("phi" = "Overdispersion",
"beta0" = "Intercept",
"beta1" = "BWQS term",
"delta1" = "Testing ratio: spline term 1",
"delta2" = "Testing ratio: spline term 2",
"delta3" = "Testing ratio: spline term 3",
labels1)
# Supplemental Table 2: Parameter estimates, credible intervals, and diagnostics from main BWQS infections model
BWQS_params %>% bind_cols(., "terms" = labels2) %>%
dplyr::select(-label) %>%
mutate_at(vars(1:6), ~round(., 3)) %>%
mutate_at(vars(5), ~round(., 0)) %>%
mutate(`95% CrI`=paste0("(",format(`2.5%`, nsmall=3),", ",format(`97.5%`, nsmall=3),")")) %>%
mutate(terms =
case_when(terms=="BWQS term" ~ "COVID-19 inequity index",
TRUE ~ terms)) %>%
dplyr::select(-`2.5%`, -`97.5%`) %>%
rename(median=`50%`) %>%
dplyr::select(terms, mean, `95% CrI`, median, Rhat, n_eff) %>%
kbl(align=rep('r', 6), font_size = 9.5) %>%
kable_classic(full_width = F, html_font = "Arial")| terms | mean | 95% CrI | median | Rhat | n_eff |
|---|---|---|---|---|---|
| Overdispersion | 104.271 | (81.793, 129.364) | 103.650 | 1.000 | 1730 |
| Intercept | 6.261 | ( 6.177, 6.344) | 6.261 | 0.999 | 1663 |
| COVID-19 inequity index | 0.074 | ( 0.063, 0.086) | 0.074 | 1.001 | 1631 |
| Testing ratio: spline term 1 | 0.983 | ( 0.893, 1.072) | 0.985 | 1.000 | 1781 |
| Testing ratio: spline term 2 | 2.029 | ( 1.821, 2.244) | 2.025 | 1.000 | 1539 |
| Testing ratio: spline term 3 | 1.121 | ( 1.002, 1.233) | 1.121 | 0.999 | 1546 |
| 1/ Median income | 0.105 | ( 0.009, 0.231) | 0.100 | 1.000 | 1638 |
| Uninsured | 0.177 | ( 0.056, 0.302) | 0.176 | 1.000 | 1505 |
| Unemployed | 0.027 | ( 0.001, 0.076) | 0.022 | 1.000 | 1672 |
| 1/ Grocers per 1000 | 0.055 | ( 0.003, 0.143) | 0.049 | 1.000 | 1660 |
| Essential Workers | 0.062 | ( 0.002, 0.173) | 0.052 | 0.999 | 1528 |
| Essential Worker: Public Transit | 0.050 | ( 0.002, 0.138) | 0.042 | 0.999 | 1771 |
| Essential Worker: Driving Commute | 0.166 | ( 0.035, 0.289) | 0.167 | 1.000 | 1790 |
| 1/ Work from home | 0.089 | ( 0.009, 0.193) | 0.086 | 1.000 | 1706 |
| Population Density by Residential Volume | 0.108 | ( 0.013, 0.210) | 0.107 | 0.999 | 1798 |
| Average people per household | 0.161 | ( 0.032, 0.297) | 0.158 | 0.999 | 1715 |
# create a figure with parameter estimates for the weights
fig2 <- ggplot(data=BWQS_fits, aes(x= reorder(label, mean), y=mean, ymin=lower, ymax=upper)) +
geom_pointrange() +
coord_flip() +
xlab("") +
ylab("Mean (95% credible intervals)") +
scale_x_discrete(labels = labels1) +
theme_set(theme_bw(base_size = 18)) +
facet_grid(group~., scales = "free", space = "free") +
theme(strip.text.x = element_text(size = 14))
if(export.figs) ggsave(fig2, filename = file.path(fig.path, paste0("fig2", "_", Sys.Date(),".png")),
dpi = 300, width = 8, height = 8)
# Step 4: Construct the summative COVID-19 inequity index value for each ZCTA
# Use the variable-specific weight on the quantiled splits to create a 10 point ZCTA-level infection risk score
BWQS_weights <- as.numeric(summary(m1)$summary[(length(vars)-length(SES_vars) + 1):number_of_coefficients,c(1)])
ZCTA_ACS_COVID2 <- m1data$X * BWQS_weights[col(ZCTA_ACS_COVID)]
BWQS_index <- ZCTA_ACS_COVID2 %>%
dplyr::mutate(BWQS_index = rowSums(.)) %>%
dplyr::select(BWQS_index)
# Visualize the relationship between BWQS and test_ratio
ggplot(data.frame(BWQS_index = BWQS_index$BWQS_index, testing_ratio = ZCTA_ACS_COVID$testing_ratio), aes(testing_ratio, BWQS_index)) + geom_point() +
geom_smooth() + geom_smooth(color = "red", method = "lm")## `geom_smooth()` using method = 'loess' and formula 'y ~ x'## `geom_smooth()` using formula 'y ~ x'
# BWQS is correlated with testing_ratio
cor.test(BWQS_index$BWQS_index, ZCTA_ACS_COVID$testing_ratio, method = "spearman") ##
## Spearman's rank correlation rho
##
## data: BWQS_index$BWQS_index and ZCTA_ACS_COVID$testing_ratio
## S = 472628, p-value = 1.969e-12
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4885952# Partial regression plot for BWQS index
# shows a nice linear relationship of BWQS index with infections after adjustment for testing
nb_testing_ns3df <- glm.nb(m1data$data_list$y ~ m1data$K)
yresid <- resid(nb_testing_ns3df)
bwqs_testing_ns3df <- lm(BWQS_index$BWQS_index ~ m1data$K)
bwqsresid <- resid(bwqs_testing_ns3df)
ggplot(data.frame(yresid, bwqsresid), aes(bwqsresid, yresid)) + geom_point() +
geom_smooth(formula = y ~ x, method = "lm", se = F) +
ylab("residual log infection rate\n(adjusted for testing)") + xlab("residual BWQS infection risk index\n(adjusted for testing)")
summary(lm(yresid ~ bwqsresid))##
## Call:
## lm(formula = yresid ~ bwqsresid)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.61379 -0.36053 0.01232 0.32793 1.95081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.05117 0.04794 -1.067 0.287
## bwqsresid 0.48714 0.03001 16.230 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6378 on 175 degrees of freedom
## Multiple R-squared: 0.6008, Adjusted R-squared: 0.5986
## F-statistic: 263.4 on 1 and 175 DF, p-value: < 2.2e-16rm(nb_testing_ns3df, yresid, bwqs_testing_ns3df, bwqsresid)
# Visualize the full model observed vs predicted -- shows a close relationship (which is why R2 is so high)
preddf <- data.frame(data.frame(pred = colMeans(extract(m1,"y_new")$y_new), y = m1data$data_list$y))
ggplot(preddf, aes(pred, y)) +
geom_point() +
geom_abline(linetype = "dashed", color = "grey10") +
coord_fixed(ratio = 1) +
scale_x_continuous("Predicted Infections per 100,000", label=comma, limits = range(preddf)) +
scale_y_continuous("Infections per 100,000", label=comma) +
theme_bw()
rm(preddf)
# calculate credible interval over the mean predicted infections
# at the median testing_ratio
sim_out <- data.frame(extract(m1, pars = c("beta0", "beta1", "delta")))
median_testing <- predict(m1data$K, median(ZCTA_ACS_COVID$testing_ratio))
# calculate term for median testing_rate
sim_out$deltamedian <- with(sim_out, median_testing[1] * delta.1 +
median_testing[2] * delta.2 +
median_testing[3] * delta.3)
# make a sequence of BWQS values
xseqlength <- 300
bwqs_seq <- seq(from = min(BWQS_index$BWQS_index),
to = max(BWQS_index$BWQS_index),
length.out = xseqlength)
sim_matrix <- sim_out$beta0 %o% rep(1, xseqlength) +
sim_out$beta1 %o% bwqs_seq +
sim_out$deltamedian%o% rep(1, xseqlength)
sim_bwqs_df <- data.frame(bwqs_seq,
lower = exp(colQuantiles(sim_matrix, probs = 0.025)),
upper = exp(colQuantiles(sim_matrix, probs = 0.975)),
mean = exp(colMeans(sim_matrix)))
rm(median_testing, xseqlength, bwqs_seq, sim_matrix)
# calculate credible interval over the mean predicted infections
# at the median COVID-19 inequity index
sim_out$beta1median <- sim_out$beta1 * median(BWQS_index$BWQS_index)
# make a sequence of testing_ratio values
xseqlength <- 300
inequity_seq <- data.frame(x = seq(from = min(ZCTA_ACS_COVID$testing_ratio),
to = max(ZCTA_ACS_COVID$testing_ratio),
length.out = xseqlength))
inequity_seq <- data.frame(x = inequity_seq$x, t(sapply(inequity_seq$x, FUN = function(x) predict(m1data$K, x))))
sim_matrix <-
sim_out$beta0 %o% rep(1, xseqlength) +
sim_out$beta1median %o% rep(1, xseqlength) +
sim_out$delta.1 %o% inequity_seq$X1 +
sim_out$delta.2 %o% inequity_seq$X2 +
sim_out$delta.3 %o% inequity_seq$X3
sim_testing_df <- data.frame(x = inequity_seq$x,
lower = exp(colQuantiles(sim_matrix, probs = 0.025)),
upper = exp(colQuantiles(sim_matrix, probs = 0.975)),
mean = exp(colMeans(sim_matrix)))
rm(sim_out, xseqlength, inequity_seq, sim_matrix)
# Visualize the relationship between BWQS index and infection rate at the median testing_ratio
BWQS_scatter <- ggplot(sim_bwqs_df, aes(x = bwqs_seq)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = 'grey75') +
geom_line(aes(y = mean)) +
geom_point(data = data.frame(BWQS_index, y = m1data$data_list$y), aes(BWQS_index, y), alpha = 0.5) +
scale_x_continuous("COVID-19 inequity index") +
scale_y_continuous("Infections per 100,000", label=comma)
BWQS_scatter <- ggExtra::ggMarginal(BWQS_scatter, type = "histogram", fill = "grey40", margins = "both",
xparams = list(binwidth = 0.5), yparams = list(binwidth = 200))
if(export.figs) {
png(filename = file.path(fig.path, paste0("fig1_", Sys.Date(), ".png")), width = 96*5, height = 96*5)
print(BWQS_scatter)
dev.off()
}## png
## 2
# Visualize the relationship between testing_ratio and infection rate at the median BWQS
testing_scatter <- ggplot(sim_testing_df, aes(x = x)) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = 'grey75') +
geom_line(aes(y = mean)) +
geom_point(data = data.frame(testing_ratio = ZCTA_ACS_COVID$testing_ratio, y = m1data$data_list$y),
aes(testing_ratio, y), alpha = 0.5) +
scale_x_continuous("Testing ratio") +
scale_y_continuous("Infections per 100,000", label=comma)
testing_scatter <- ggExtra::ggMarginal(testing_scatter, type = "histogram", fill = "grey40", margins = "both",
xparams = list(binwidth = 0.005), yparams = list(binwidth = 200))
if(export.figs) {
png(filename = file.path(fig.path, paste0("sfig2_", Sys.Date(), ".png")), width = 96*5, height = 96*5)
print(testing_scatter)
dev.off()
}## png
## 2
# Step 5: Visualize the spatial distribution of ZCTA-level infection risk scores
ZCTA_BWQS_COVID_shp <- ZCTA_ACS_COVID_shp %>%
bind_cols(., BWQS_index)
# reproject to WGS84 to be compatible with scalebar
ZCTA_BWQS_COVID_shp <- st_transform(ZCTA_BWQS_COVID_shp, 4326)
# Figure 3A - BWQS Index by ZCTA
fig3a <- ggplot() +
geom_sf(data = basemap_water, fill = "white", lwd = 0) +
geom_sf(data = ZCTA_BWQS_COVID_shp, aes(fill = BWQS_index), lwd = 0.2) +
scalebar(ZCTA_BWQS_COVID_shp, dist = 5, dist_unit = "km",
transform = TRUE, model = "WGS84",
st.size = 2.8, height = 0.015, border.size = 0.5, st.dist = 0.011,
anchor = c(x = -73.71, y = 40.49)) +
scale_fill_gradientn(colours=brewer_pal("YlGnBu", type = "seq")(7)) +
coord_sf(crs = st_crs(ZCTA_BWQS_COVID_shp),
xlim = c(plot_bounds$xmin, plot_bounds$xmax),
ylim = c(plot_bounds$ymin, plot_bounds$ymax),
expand = FALSE) +
theme_bw(base_size = 6) +
labs(fill = "COVID-19 Inequity Index") +
theme(legend.title = element_text(face = "bold", size = 9),
panel.background = element_rect(fill = "#cccccc"),
legend.background = element_rect(fill = "transparent"),
panel.grid = element_blank(),
legend.position = c(0.125, 0.90),
legend.text = element_text(size = 8.5),
legend.key.size = unit(1.1, "lines"),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_blank(),
axis.title.x = element_blank())
# Step 6: Compare quantile distribution of ZCTA-level BWQS scores by the race/ethnic composition of residents
Demographics <- ACS_Data_scaled %>%
dplyr::select(GEOID, ends_with("_raceethnic"), total_pop1) %>%
left_join(., modzcta_to_zcta1, by = c("GEOID" = "ZCTA")) %>%
group_by(MODZCTA) %>%
summarise_at(vars(ends_with("_raceethnic"), total_pop1), ~sum(.)) %>%
mutate(zcta = as.character(MODZCTA),
other_raceethnic = total_pop1 - (hisplat_raceethnic + nonhispLat_white_raceethnic + nonhispLat_black_raceethnic +
nonhispLat_amerindian_raceethnic + nonhispLat_asian_raceethnic))
ZCTA_BQWS <- ZCTA_BWQS_COVID_shp %>%
st_set_geometry(., NULL) %>%
dplyr::select(zcta, BWQS_index)
Demographics_for_ridges <- Demographics %>%
left_join(., ZCTA_BQWS, by = "zcta") %>%
dplyr::select(-total_pop1) %>%
gather(key = "Race/Ethnicity", value = "Population", hisplat_raceethnic:other_raceethnic) %>%
mutate(`Race/Ethnicity` = if_else(`Race/Ethnicity`=="hisplat_raceethnic","Hispanic/Latinx",
if_else(`Race/Ethnicity`=="nonhispLat_black_raceethnic", "Black",
if_else(`Race/Ethnicity`=="nonhispLat_white_raceethnic", "White",
if_else(`Race/Ethnicity`== "nonhispLat_asian_raceethnic", "Asian", "Other")))),
`Race/Ethnicity` = factor(`Race/Ethnicity`, levels = c( "White", "Asian", "Other","Hispanic/Latinx","Black")),
Population = as.numeric(Population))
(Demographics_for_ridges %>%
group_by(`Race/Ethnicity`) %>%
summarise(weighted.mean(BWQS_index, Population),
weightedMedian(BWQS_index, Population)))## # A tibble: 5 x 3
## `Race/Ethnicity` `weighted.mean(BWQS_index, Po… `weightedMedian(BWQS_index, P…
## <fct> <dbl> <dbl>
## 1 White 4.34 4.71
## 2 Asian 5.48 6.20
## 3 Other 5.14 5.67
## 4 Hispanic/Latinx 6.10 6.60
## 5 Black 6.26 6.51fig4 <- ggplot(Demographics_for_ridges,
aes(x = BWQS_index, y = `Race/Ethnicity`)) +
xlab("COVID-19 inequity index")+
geom_density_ridges(
aes(height = ..density..,
weight = Population / sum(Population)),
scale = 0.95,
stat ="density") +
scale_y_discrete(expand = c(0, 0)) +
expand_limits(y = 6) +
theme(legend.position = "none")
if(export.figs) ggsave(plot = fig4, filename = file.path(fig.path, paste0("fig4","_",Sys.Date(),".png")),
dpi = 300, device = "png", width = 8, height = 5)
Below_25th_zctas <- ZCTA_BQWS %>%
filter(BWQS_index<quantile(BWQS_index, .25))Race_Ethncity_below25th <- Demographics %>%
filter(zcta %in% Below_25th_zctas$zcta) %>%
dplyr::select(-zcta, -MODZCTA) %>%
summarise_all(funs(sum(.))) %>%
mutate_at(vars(ends_with("_raceethnic")), funs(round((./total_pop1)*100, 2))) %>%
mutate(Group = "Below 25th percentile BWQS")
Between_25thand75th_zctas <- ZCTA_BQWS %>%
filter(BWQS_index>quantile(BWQS_index, .25) & BWQS_index<quantile(BWQS_index, .75))
Race_Ethncity_btw_25th75th <- Demographics %>%
filter(zcta %in% Between_25thand75th_zctas$zcta) %>%
dplyr::select(-zcta, -MODZCTA) %>%
summarise_all(funs(sum(.))) %>%
mutate_at(vars(ends_with("_raceethnic")), funs(round((./total_pop1)*100, 2)))%>%
mutate(Group = "Between 25-75th percentile BWQS")
Above_75th_zctas <- ZCTA_BQWS %>%
filter(BWQS_index>quantile(BWQS_index, .75))
Race_Ethncity_above_75th <- Demographics %>%
filter(zcta %in% Above_75th_zctas$zcta) %>%
dplyr::select(-zcta, -MODZCTA) %>%
summarise_all(funs(sum(.))) %>%
mutate_at(vars(ends_with("_raceethnic")), funs(round((./total_pop1)*100, 2))) %>%
mutate(Group = "Above 75th percentile BWQS")
Race_Ethncity_all <- Demographics %>%
dplyr::select(-zcta, -MODZCTA) %>%
summarise_all(funs(sum(.))) %>%
mutate_at(vars(ends_with("_raceethnic")), funs(round((./total_pop1)*100, 2)))%>%
mutate(Group = "NYC Population")
Demographics_by_BWQS <- bind_rows(Race_Ethncity_all, Race_Ethncity_below25th, Race_Ethncity_btw_25th75th, Race_Ethncity_above_75th) %>%
mutate(Other = other_raceethnic + nonhispLat_amerindian_raceethnic) %>%
dplyr::select(-other_raceethnic, -nonhispLat_amerindian_raceethnic, - total_pop1) %>%
rename("Hispanic/Latinx" = "hisplat_raceethnic",
"Black" = "nonhispLat_black_raceethnic",
"White" = "nonhispLat_white_raceethnic",
"Asian" = "nonhispLat_asian_raceethnic") %>%
gather(., "Race/Ethnicity", "Proportion", 1:4,6) %>%
mutate(`Race/Ethnicity` = factor(`Race/Ethnicity`, levels = c("Other", "Asian", "White", "Hispanic/Latinx", "Black")),
Group = factor(Group, levels = c("NYC Population", "Below 25th percentile BWQS", "Between 25-75th percentile BWQS", "Above 75th percentile BWQS")))
labels_demographics <- c("NYC Population" = "NYC Population", "Below 25th percentile BWQS" = "Below 25th\npercentile BWQS",
"Between 25-75th percentile BWQS" = "Between 25-75th\npercentile BWQS", "Above 75th percentile BWQS" = "Above 75th\npercentile BWQS")
labels_demographics <- c("NYC Population" = "NYC Population", "Below 25th percentile BWQS" = "Below 25th\n%ile index",
"Between 25-75th percentile BWQS" = "Between 25-75th\n%ile index", "Above 75th percentile BWQS" = "Above 75th\n%ile index")
sfig4 <- ggplot(Demographics_by_BWQS, aes(fill=`Race/Ethnicity`, y=Proportion, x=Group)) +
geom_rect(data = subset(Demographics_by_BWQS, Group=="NYC Population"),
aes(xmin=as.numeric(Group)-.35,xmax=as.numeric(Group)+.35, ymin=0, ymax=100, fill="gray85"), color = "gray", alpha = .1) +
geom_bar(data = subset(Demographics_by_BWQS, Group!="NYC Population"), position="stack", stat="identity", width = .75) +
geom_bar(data = subset(Demographics_by_BWQS, Group=="NYC Population"), position="stack", stat="identity", width = .45) +
scale_fill_manual(breaks = c("Other", "Asian", "White", "Hispanic/Latinx", "Black"),
values = c("#984ea3","#ff7f00","gray85", "#4daf4a", "#e41a1c", "#377eb8")) +
geom_text(aes(label=ifelse(Proportion >= 5, paste0(sprintf("%.0f", Proportion),"%"),"")),
position=position_stack(vjust=0.5), colour="black", size = 8) +
scale_y_continuous("Percentage", expand = c(0,0), labels = function(x) paste0(x, "%")) +
scale_x_discrete(limits = c( "NYC Population", "Below 25th percentile BWQS", "Between 25-75th percentile BWQS", "Above 75th percentile BWQS"),
labels = labels_demographics) +
theme_bw(base_size = 16) +
theme(legend.text = element_text(size = 14),
axis.text.y = element_text(size = 16),
axis.text.x = element_text(size = 16, color = "black"),
axis.title.x = element_blank())
if(export.figs) ggsave(sfig4, filename = file.path(fig.path, paste0("sfig4","_",Sys.Date(),".png")), device = "png",
dpi = 300, width = 12, height = 6)
t_postm1 = Sys.time()BWQS by Tract
#### BWQS by Tract ####
get_tract_vars_by_zcta <- function(tract_vars, modzcta_tract_crosswalk, whichq){
whichq = str_to_lower(whichq)
if(whichq %in% c("median", "q2", "2q", "2")){
qname = "median"
qnum = 2
} else if(whichq %in% c("q3", "3q", "3")){
qname = "3Q"
qnum = 3
} else { stop("unhandled quartile selected: ", whichq)}
ses_zcta <- tract_vars %>%
st_set_geometry(., NULL) %>%
left_join(., modzcta_tract_crosswalk, by = c("GEOID" = "TRACT")) %>%
filter(!is.na(MODZCTA)) %>%
group_by(MODZCTA) %>%
summarise(essentialworker_drove_rename = Hmisc::wtd.quantile(essentialworker_drove, W_RES_RATIO)[[qnum]],
essentialworker_pubtrans_rename = Hmisc::wtd.quantile(essentialworker_pubtrans, W_RES_RATIO)[[qnum]],
not_quarantined_jobs_rename = Hmisc::wtd.quantile(not_quarantined_jobs, W_RES_RATIO)[[qnum]],
didnot_workhome_commute_rename = Hmisc::wtd.quantile(didnot_workhome_commute, W_RES_RATIO)[[qnum]],
not_insured_rename = Hmisc::wtd.quantile(not_insured, W_RES_RATIO)[[qnum]],
one_over_medincome_rename = Hmisc::wtd.quantile(one_over_medincome, W_RES_RATIO)[[qnum]],
unemployed_rename = Hmisc::wtd.quantile(unemployed, W_RES_RATIO)[[qnum]],
avg_hhold_size_rename = Hmisc::wtd.quantile(avg_hhold_size, W_RES_RATIO)[[qnum]],
res_vol_popdensity_rename = Hmisc::wtd.quantile(res_vol_popdensity, W_RES_RATIO)[[qnum]],
one_over_grocers_per_1000_rename = Hmisc::wtd.quantile(one_over_grocers_per_1000, W_RES_RATIO)[[qnum]])
ses_zcta %>% rename_with(~ gsub("rename", qname, .x))
}
# Select median tract values by MODZCTA, weighted by residential address share and population density by tract
SES_zcta_median <- get_tract_vars_by_zcta(tract_vars, modzcta_to_tract2, "median")
SES_vars_median = names(SES_zcta_median)[names(SES_zcta_median) != "MODZCTA"]
# join SES to testing data: positive/100k and testing_ratio
SES_zcta_median_testing <- ZCTA_ACS_COVID %>%
dplyr::select(zcta, pos_per_100000, testing_ratio) %>%
left_join(SES_zcta_median, by = c("zcta" = "MODZCTA"))
# BWQS using median at the tract level
m2data_median <- prep_BWQS_data(SES_zcta_median_testing, SES_vars_median)
m2_median <- fit_BWQS_model(data_list = m2data_median$data_list, stan_model_path = BWQS_stan_model)## Warning in validityMethod(object): The following variables have undefined
## values: log_lik[1],The following variables have undefined values: log_lik[2],The
## following variables have undefined values: log_lik[3],The following variables
## have undefined values: log_lik[4],The following variables have undefined
## values: log_lik[5],The following variables have undefined values: log_lik[6],The
## following variables have undefined values: log_lik[7],The following variables
## have undefined values: log_lik[8],The following variables have undefined values:
## log_lik[9],The following variables have undefined values: log_lik[10],The
## following variables have undefined values: log_lik[11],The following
## variables have undefined values: log_lik[12],The following variables have
## undefined values: log_lik[13],The following variables have undefined values:
## log_lik[14],The following variables have undefined values: log_lik[15],The
## following variables have undefined values: log_lik[16],The following
## variables have undefined values: log_lik[17],The following variables have
## undefined values: log_lik[18],The following variables have undefined values:
## log_lik[19],The following variables have undefined values: log_lik[20],The
## following variables have undefined values: log_lik[21],The following
## variables have undefined values: log_lik[22],The following variables have
## undefined values: log_lik[23],The following variables have undefined values:
## log_lik[24],The following variables have undefined values: log_lik[25],The
## following variables have undefined values: log_lik[26],The following
## variables have undefined values: log_lik[27],The following variables have
## undefined values: log_lik[28],The following variables have undefined values:
## log_lik[29],The following variables have undefined values: log_lik[30],The
## following variables have undefined values: log_lik[31],The following
## variables have undefined values: log_lik[32],The following variables have
## undefined values: log_lik[33],The following variables have undefined values:
## log_lik[34],The following variables have undefined values: log_lik[35],The
## following variables have undefined values: log_lik[36],The following
## variables have undefined values: log_lik[37],The following variables have
## undefined values: log_lik[38],The following variables have undefined values:
## log_lik[39],The following variables have undefined values: log_lik[40],The
## following variables have undefined values: log_lik[41],The following
## variables have undefined values: log_lik[42],The following variables have
## undefined values: log_lik[43],The following variables have undefined values:
## log_lik[44],The following variables have undefined values: log_lik[45],The
## following variables have undefined values: log_lik[46],The following
## variables have undefined values: log_lik[47],The following variables have
## undefined values: log_lik[48],The following variables have undefined values:
## log_lik[49],The following variables have undefined values: log_lik[50],The
## following variables have undefined values: log_lik[51],The following
## variables have undefined values: log_lik[52],The following variables have
## undefined values: log_lik[53],The following variables have undefined values:
## log_lik[54],The following variables have undefined values: log_lik[55],The
## following variables have undefined values: log_lik[56],The following
## variables have undefined values: log_lik[57],The following variables have
## undefined values: log_lik[58],The following variables have undefined values:
## log_lik[59],The following variables have undefined values: log_lik[60],The
## following variables have undefined values: log_lik[61],The following
## variables have undefined values: log_lik[62],The following variables have
## undefined values: log_lik[63],The following variables have undefined values:
## log_lik[64],The following variables have undefined values: log_lik[65],The
## following variables have undefined values: log_lik[66],The following
## variables have undefined values: log_lik[67],The following variables have
## undefined values: log_lik[68],The following variables have undefined values:
## log_lik[69],The following variables have undefined values: log_lik[70],The
## following variables have undefined values: log_lik[71],The following
## variables have undefined values: log_lik[72],The following variables have
## undefined values: log_lik[73],The following variables have undefined values:
## log_lik[74],The following variables have undefined values: log_lik[75],The
## following variables have undefined values: log_lik[76],The following
## variables have undefined values: log_lik[77],The following variables have
## undefined values: log_lik[78],The following variables have undefined values:
## log_lik[79],The following variables have undefined values: log_lik[80],The
## following variables have undefined values: log_lik[81],The following
## variables have undefined values: log_lik[82],The following variables have
## undefined values: log_lik[83],The following variables have undefined values:
## log_lik[84],The following variables have undefined values: log_lik[85],The
## following variables have undefined values: log_lik[86],The following
## variables have undefined values: log_lik[87],The following variables have
## undefined values: log_lik[88],The following variables have undefined values:
## log_lik[89],The following variables have undefined values: log_lik[90],The
## following variables have undefined values: log_lik[91],The following
## variables have undefined values: log_lik[92],The following variables have
## undefined values: log_lik[93],The following variables have undefined values:
## log_lik[94],The following variables have undefined values: log_lik[95],The
## following variables have undefined values: log_lik[96],The following
## variables have undefined values: log_lik[97],The following variables have
## undefined values: log_lik[98],The following variables have undefined values:
## log_lik[99],The following variables have undefined values: log_lik[100],The
## following variables have undefined values: log_lik[101],The following
## variables have undefined values: log_lik[102],The following variables have
## undefined values: log_lik[103],The following variables have undefined values:
## log_lik[104],The following variables have undefined values: log_lik[105],The
## following variables have undefined values: log_lik[106],The following
## variables have undefined values: log_lik[107],The following variables have
## undefined values: log_lik[108],The following variables have undefined values:
## log_lik[109],The following variables have undefined values: log_lik[110],The
## following variables have undefined values: log_lik[111],The following
## variables have undefined values: log_lik[112],The following variables have
## undefined values: log_lik[113],The following variables have undefined values:
## log_lik[114],The following variables have undefined values: log_lik[115],The
## following variables have undefined values: log_lik[116],The following
## variables have undefined values: log_lik[117],The following variables have
## undefined values: log_lik[118],The following variables have undefined values:
## log_lik[119],The following variables have undefined values: log_lik[120],The
## following variables have undefined values: log_lik[121],The following
## variables have undefined values: log_lik[122],The following variables have
## undefined values: log_lik[123],The following variables have undefined values:
## log_lik[124],The following variables have undefined values: log_lik[125],The
## following variables have undefined values: log_lik[126],The following
## variables have undefined values: log_lik[127],The following variables have
## undefined values: log_lik[128],The following variables have undefined values:
## log_lik[129],The following variables have undefined values: log_lik[130],The
## following variables have undefined values: log_lik[131],The following
## variables have undefined values: log_lik[132],The following variables have
## undefined values: log_lik[133],The following variables have undefined values:
## log_lik[134],The following variables have undefined values: log_lik[135],The
## following variables have undefined values: log_lik[136],Tht_m2 = Sys.time()
# Select third quartile tract values by MODZCTA, weighted by residential address share and population density by tract
SES_zcta_3q <- get_tract_vars_by_zcta(tract_vars, modzcta_to_tract2, "3q")
SES_vars_3q = names(SES_zcta_3q)[names(SES_zcta_3q) != "MODZCTA"]
SES_zcta_3q_testing <- ZCTA_ACS_COVID %>%
dplyr::select(zcta, pos_per_100000, testing_ratio) %>%
left_join(SES_zcta_3q, by = c("zcta" = "MODZCTA"))
#BWQS using the 3rd quartile at the tract level
m2data_3q <- prep_BWQS_data(SES_zcta_3q_testing, SES_vars_3q)
m2_3q <- fit_BWQS_model(data_list = m2data_3q$data_list, stan_model_path = BWQS_stan_model)## Warning in validityMethod(object): The following variables have undefined
## values: log_lik[1],The following variables have undefined values: log_lik[2],The
## following variables have undefined values: log_lik[3],The following variables
## have undefined values: log_lik[4],The following variables have undefined
## values: log_lik[5],The following variables have undefined values: log_lik[6],The
## following variables have undefined values: log_lik[7],The following variables
## have undefined values: log_lik[8],The following variables have undefined values:
## log_lik[9],The following variables have undefined values: log_lik[10],The
## following variables have undefined values: log_lik[11],The following
## variables have undefined values: log_lik[12],The following variables have
## undefined values: log_lik[13],The following variables have undefined values:
## log_lik[14],The following variables have undefined values: log_lik[15],The
## following variables have undefined values: log_lik[16],The following
## variables have undefined values: log_lik[17],The following variables have
## undefined values: log_lik[18],The following variables have undefined values:
## log_lik[19],The following variables have undefined values: log_lik[20],The
## following variables have undefined values: log_lik[21],The following
## variables have undefined values: log_lik[22],The following variables have
## undefined values: log_lik[23],The following variables have undefined values:
## log_lik[24],The following variables have undefined values: log_lik[25],The
## following variables have undefined values: log_lik[26],The following
## variables have undefined values: log_lik[27],The following variables have
## undefined values: log_lik[28],The following variables have undefined values:
## log_lik[29],The following variables have undefined values: log_lik[30],The
## following variables have undefined values: log_lik[31],The following
## variables have undefined values: log_lik[32],The following variables have
## undefined values: log_lik[33],The following variables have undefined values:
## log_lik[34],The following variables have undefined values: log_lik[35],The
## following variables have undefined values: log_lik[36],The following
## variables have undefined values: log_lik[37],The following variables have
## undefined values: log_lik[38],The following variables have undefined values:
## log_lik[39],The following variables have undefined values: log_lik[40],The
## following variables have undefined values: log_lik[41],The following
## variables have undefined values: log_lik[42],The following variables have
## undefined values: log_lik[43],The following variables have undefined values:
## log_lik[44],The following variables have undefined values: log_lik[45],The
## following variables have undefined values: log_lik[46],The following
## variables have undefined values: log_lik[47],The following variables have
## undefined values: log_lik[48],The following variables have undefined values:
## log_lik[49],The following variables have undefined values: log_lik[50],The
## following variables have undefined values: log_lik[51],The following
## variables have undefined values: log_lik[52],The following variables have
## undefined values: log_lik[53],The following variables have undefined values:
## log_lik[54],The following variables have undefined values: log_lik[55],The
## following variables have undefined values: log_lik[56],The following
## variables have undefined values: log_lik[57],The following variables have
## undefined values: log_lik[58],The following variables have undefined values:
## log_lik[59],The following variables have undefined values: log_lik[60],The
## following variables have undefined values: log_lik[61],The following
## variables have undefined values: log_lik[62],The following variables have
## undefined values: log_lik[63],The following variables have undefined values:
## log_lik[64],The following variables have undefined values: log_lik[65],The
## following variables have undefined values: log_lik[66],The following
## variables have undefined values: log_lik[67],The following variables have
## undefined values: log_lik[68],The following variables have undefined values:
## log_lik[69],The following variables have undefined values: log_lik[70],The
## following variables have undefined values: log_lik[71],The following
## variables have undefined values: log_lik[72],The following variables have
## undefined values: log_lik[73],The following variables have undefined values:
## log_lik[74],The following variables have undefined values: log_lik[75],The
## following variables have undefined values: log_lik[76],The following
## variables have undefined values: log_lik[77],The following variables have
## undefined values: log_lik[78],The following variables have undefined values:
## log_lik[79],The following variables have undefined values: log_lik[80],The
## following variables have undefined values: log_lik[81],The following
## variables have undefined values: log_lik[82],The following variables have
## undefined values: log_lik[83],The following variables have undefined values:
## log_lik[84],The following variables have undefined values: log_lik[85],The
## following variables have undefined values: log_lik[86],The following
## variables have undefined values: log_lik[87],The following variables have
## undefined values: log_lik[88],The following variables have undefined values:
## log_lik[89],The following variables have undefined values: log_lik[90],The
## following variables have undefined values: log_lik[91],The following
## variables have undefined values: log_lik[92],The following variables have
## undefined values: log_lik[93],The following variables have undefined values:
## log_lik[94],The following variables have undefined values: log_lik[95],The
## following variables have undefined values: log_lik[96],The following
## variables have undefined values: log_lik[97],The following variables have
## undefined values: log_lik[98],The following variables have undefined values:
## log_lik[99],The following variables have undefined values: log_lik[100],The
## following variables have undefined values: log_lik[101],The following
## variables have undefined values: log_lik[102],The following variables have
## undefined values: log_lik[103],The following variables have undefined values:
## log_lik[104],The following variables have undefined values: log_lik[105],The
## following variables have undefined values: log_lik[106],The following
## variables have undefined values: log_lik[107],The following variables have
## undefined values: log_lik[108],The following variables have undefined values:
## log_lik[109],The following variables have undefined values: log_lik[110],The
## following variables have undefined values: log_lik[111],The following
## variables have undefined values: log_lik[112],The following variables have
## undefined values: log_lik[113],The following variables have undefined values:
## log_lik[114],The following variables have undefined values: log_lik[115],The
## following variables have undefined values: log_lik[116],The following
## variables have undefined values: log_lik[117],The following variables have
## undefined values: log_lik[118],The following variables have undefined values:
## log_lik[119],The following variables have undefined values: log_lik[120],The
## following variables have undefined values: log_lik[121],The following
## variables have undefined values: log_lik[122],The following variables have
## undefined values: log_lik[123],The following variables have undefined values:
## log_lik[124],The following variables have undefined values: log_lik[125],The
## following variables have undefined values: log_lik[126],The following
## variables have undefined values: log_lik[127],The following variables have
## undefined values: log_lik[128],The following variables have undefined values:
## log_lik[129],The following variables have undefined values: log_lik[130],The
## following variables have undefined values: log_lik[131],The following
## variables have undefined values: log_lik[132],The following variables have
## undefined values: log_lik[133],The following variables have undefined values:
## log_lik[134],The following variables have undefined values: log_lik[135],The
## following variables have undefined values: log_lik[136],Tht_m3 = Sys.time()
Summarize_BWQS_performance <- function(model, df){
model_type = deparse(substitute(model))
waic = extract_waic(model)$waic[[1]]
bayesr2 = Compute_Bayes_R2(model)$meanr2
rmse = sqrt(mean((df$pos_per_100000 - colMeans(extract(model,"y_new")$y_new))^2))
effect_estimate = exp(mean(extract(model, "beta1")$beta1))
summarized <- bind_cols(model_name = as.character(model_type), waic = waic, bayesr2 = bayesr2, rmse = rmse, effect_estimate = effect_estimate)
return(summarized)
}
bind_rows(Summarize_BWQS_performance(m1, SES_zcta_median_testing),
Summarize_BWQS_performance(m2_median, SES_zcta_median_testing),
Summarize_BWQS_performance(m2_3q, SES_zcta_3q_testing),
)## # A tibble: 3 x 5
## model_name waic bayesr2 rmse effect_estimate
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 m1 2370. 0.934 184. 1.08
## 2 m2_median 2362. 0.936 185. 1.08
## 3 m2_3q 2359. 0.937 183. 1.08Part 2: Compare capacity to socially distance (as measured by transit data) by neighborhood-level risk scores
#### Part 2: Compare capacity to socially distance (as measured by transit data) by neighborhood-level risk scores ####
UHF_ZipCodes1 <- UHF_ZipCodes %>%
mutate(zcta = as.character(zip),
uhf = as.character(code)) %>%
dplyr::select(zcta, uhf)
ZCTA_BWQS_score <- ZCTA_BWQS_COVID_shp %>%
st_drop_geometry() %>%
dplyr::select(zcta, BWQS_index)
UHF_BWQS <- ZCTA_ACS_COVID %>%
left_join(., UHF_ZipCodes1, by = "zcta") %>%
left_join(., ZCTA_BWQS_score, join = "zcta") %>%
group_by(uhf) %>%
summarise(uhf_weighted_bwqs = weighted.mean(BWQS_index, total_pop1)) #population weighting## Joining, by = "zcta"UHF_BWQS_COVID_shp <- UHF_shp %>%
mutate(uhf = as.character(UHFCODE)) %>%
left_join(., UHF_BWQS, by = "uhf") %>%
mutate(Risk = factor(ifelse(uhf_weighted_bwqs > median(uhf_weighted_bwqs, na.rm = T), #median split of neighborhood risk
"High", "Low"), levels = c("High", "Low")))
ggplot(UHF_BWQS_COVID_shp) +
geom_sf(aes(fill = uhf_weighted_bwqs))
Subway_ridership_by_UHF %>% filter(place=="all" & date >= "2020-01-29" & date <= "2020-04-30") %>%
ggplot() + geom_point(aes(date, usage.median.ratio, color = place))
Mean_Ridership <- Subway_ridership_by_UHF %>% #Use the mean ridership to identify the proper function for a nonlinear model
filter(date >= "2020-01-29" & date <= "2020-04-30" & place=="all") %>%
arrange(date) %>%
mutate(time_index = time(date))
Mean_Ridership %>%
ggplot() + geom_point(aes(date, usage.median.ratio))
#The weibull probability distribution works best for these data
fit_drm_w2.4 <- drm(usage.median.ratio ~ time_index, fct = W2.4(), data = Mean_Ridership)
# suppress warning about vector recycling in predict.drc.R
handler <- function(w) if( any( grepl( "Recycling array of length 1 in array-vector arithmetic is deprecated.", w) ) ){
invokeRestart( "muffleWarning" ) }
DRM_mean_predictions <- bind_cols(Mean_Ridership,
as_tibble(withCallingHandlers(predict(fit_drm_w2.4, interval = "confidence"), warning = handler )))
# Supplementary Figure 6
sfig6 <- ggplot() + geom_point(data = DRM_mean_predictions, aes(x = Mean_Ridership$date, y = Mean_Ridership$usage.median.ratio)) +
geom_ribbon(data = DRM_mean_predictions, aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50", alpha = .5) +
geom_line(aes(x = DRM_mean_predictions$date, y = DRM_mean_predictions$Prediction), color = "red") +
theme_bw(base_size = 16) +
xlab("Date") +
scale_y_continuous("Relative Subway Ridership (%)", labels = scales::percent)
if(export.figs) ggsave(sfig6, filename = file.path(fig.path, paste0("sfig6", "_", Sys.Date(), ".png")),
device = "png", dpi = 300, width = 8, height = 5)
# create a dataframe for the analysis
service_changes_in_lowsubway_areas <- tibble(date = as.Date(c("2020-02-01", "2020-02-02", "2020-02-08", "2020-02-09", "2020-02-15", "2020-02-16", "2020-02-22", "2020-02-23", "2020-02-29", "2020-03-01", "2020-03-07", "2020-03-08",
"2020-02-01", "2020-02-02", "2020-02-08", "2020-02-09", "2020-02-16", "2020-02-16")),
place = c("403","403","403","403","403","403","403","403","403","403","403","403",
"502","502","502","502","502","502"))
Subway_BWQS_df <- Subway_ridership_by_UHF %>%
left_join(., st_set_geometry(UHF_BWQS_COVID_shp, NULL), by = c("place" = "uhf")) %>%
filter(!is.na(place) & date >= "2020-01-29" & date <= "2020-04-30") %>%
group_by(place) %>%
arrange(date) %>%
mutate(time_index = time(date)) %>%
filter(!is.na(Risk) & date != "2020-02-17") %>% # removing Presidents' Day national holiday
anti_join(., service_changes_in_lowsubway_areas, by = c("date", "place")) %>% # removing low outliers due to service changes in low subway density areas
ungroup()
fit_drm_all <- drm(usage.median.ratio ~ time_index, fct = W2.4(), data = Subway_BWQS_df)
plot(fit_drm_all)
fit_drm_interact <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_df, separate = T) #using this for inference
fit_drm_interact1 <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_df, separate = F) #using this for plotting -- can't get separate = T to plot correctly!
fit_drm_risk_high <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_df, subset = Risk=="High")
fit_drm_risk_low <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_df, subset = Risk=="Low")
anova(fit_drm_all, fit_drm_interact) # comparing the mean only model to the interaction model ##
## 1st model
## fct: W2.4()
## 2nd model
## fct: NULL## ANOVA table
##
## ModelDf RSS Df F value p value
## 1st model 3291 21.3020
## 2nd model 3287 9.7062 4 981.73 0.00summary(fit_drm_interact)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:Low -11.3031622 0.3539527 -31.934 < 2.2e-16 ***
## c:Low 0.0948445 0.0032081 29.564 < 2.2e-16 ***
## d:Low 0.9771921 0.0028358 344.594 < 2.2e-16 ***
## e:Low 44.8234259 0.0974926 459.762 < 2.2e-16 ***
## b:High -10.9375092 0.3344433 -32.704 < 2.2e-16 ***
## c:High 0.1600641 0.0033505 47.773 < 2.2e-16 ***
## d:High 0.9643656 0.0027745 347.587 < 2.2e-16 ***
## e:High 46.5688477 0.1056759 440.676 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.07299479 (3287 degrees of freedom)confint(fit_drm_interact)## 2.5 % 97.5 %
## b:Low -11.99715225 -10.6091722
## c:Low 0.08855447 0.1011346
## d:Low 0.97163205 0.9827522
## e:Low 44.63227355 45.0145782
## b:High -11.59324749 -10.2817710
## c:High 0.15349477 0.1666334
## d:High 0.95892577 0.9698054
## e:High 46.36165053 46.7760449compParm(fit_drm_interact, "b", "-")##
## Comparison of parameter 'b'
##
## Estimate Std. Error t-value p-value
## Low-High -0.36565 0.48696 -0.7509 0.4528compParm(fit_drm_interact, "c", "-")##
## Comparison of parameter 'c'
##
## Estimate Std. Error t-value p-value
## Low-High -0.0652196 0.0046387 -14.06 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1slopes <- as_tibble(coef(fit_drm_interact), rownames = "vars") %>%
separate(col = vars, into = c("vars", "BWQS_risk"), sep = ":") %>%
spread(., vars, value) %>%
mutate(slope = -b*(-d/(4*e)))
as_tibble(confint(fit_drm_interact), rownames = "vars") %>%
separate(col = vars, into = c("vars", "BWQS_risk"), sep = ":") %>%
filter(vars == "b") %>%
right_join(., slopes, by = "BWQS_risk") %>%
mutate(slope_lower_ci = -`2.5 %`*(-d/(4*e)),
slope_higher_ci = -`97.5 %`*(-d/(4*e))) %>%
dplyr::select(BWQS_risk, slope, slope_lower_ci, slope_higher_ci)## # A tibble: 2 x 4
## BWQS_risk slope slope_lower_ci slope_higher_ci
## <chr> <dbl> <dbl> <dbl>
## 1 Low -0.0616 -0.0654 -0.0578
## 2 High -0.0566 -0.0600 -0.0532fit_drm_predictions <- as_tibble(withCallingHandlers(predict(fit_drm_interact1, as.data.frame(Subway_BWQS_df %>% dplyr::select(usage.median.ratio, time_index, Risk)),
interval = "confidence"), warning = handler))
Subway_BWQS_df1 <- bind_cols(Subway_BWQS_df, fit_drm_predictions)
Subway_BWQS_df2 <- Subway_BWQS_df1 %>%
filter(date>"2020-02-16") %>% # subsetting for visualization
mutate(Risk = if_else(Risk == "High", "High (above median)", "Low (below median)"))
fig5 <- ggplot() +
geom_jitter(data = Subway_BWQS_df2, aes(x = date, y = usage.median.ratio, color = Risk), alpha = .5, position = position_jitter(height = 0, width = 0.4))+
geom_ribbon(data = subset(Subway_BWQS_df2, Risk == "High (above median)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_ribbon(data = subset(Subway_BWQS_df2, Risk == "Low (below median)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_line(data = Subway_BWQS_df2, aes(x = date, y = Prediction, color = Risk)) +
scale_x_date("Date", date_minor_breaks = "1 week") +
scale_y_continuous("Relative Subway Ridership (%)", labels = scales::percent) +
geom_vline(xintercept = date("2020-03-22"),color = "grey30", lty = 2) +
theme_bw(base_size = 16) +
labs(colour="COVID-19 inequity index") +
theme(legend.title = element_text(face = "bold", size = 12), legend.position = c(0.8, 0.7))
if(export.figs) ggsave(fig5, filename = file.path(fig.path, paste0("fig5", "_", Sys.Date() ,".png")),
dpi = 300, width = 8, height = 6)
# which UHF neighborhoods were dropped?
included_uhf <- Subway_BWQS_df %>% distinct(UHFCODE)
notincluded_uhf_shp <- UHF_BWQS_COVID_shp %>%
filter(!UHFCODE %in% included_uhf$UHFCODE &
UHFCODE !=0) %>%
mutate(NotIncluded = "*")
# Supplementary Figure 5
sfig5 <- ggplot() +
geom_sf(data = basemap_water, fill = "white", lwd = 0) +
geom_sf(data = NYC_basemap_shp, size = 0.2) +
geom_sf(data = subset(UHF_BWQS_COVID_shp, !is.na(Risk)), aes(fill = Risk), size = 0.2) +
labs(fill = "COVID-19 inequity index\ngrouping") +
geom_sf(data = SubwayStation_shp, size = 0.6) +
geom_sf_text(data = notincluded_uhf_shp, aes(label = NotIncluded), size = 9, color = "grey15") +
coord_sf(crs = st_crs(ZCTA_BWQS_COVID_shp),
xlim = c(plot_bounds$xmin, plot_bounds$xmax),
ylim = c(plot_bounds$ymin, plot_bounds$ymax),
expand = FALSE) +
theme_bw() +
theme_smallmaps +
theme(legend.position = c(0.22, 0.86), plot.margin = margin(0, 0, 0, 0, "pt"))if(export.figs) ggsave(sfig5, filename = file.path(fig.path, paste0("sfig5", "_", Sys.Date(),".png")),
dpi = 300, width = 4, height = 4)
Part 3: Spatial analysis of mortality in relation to BWQS scores
#### Part 3: Spatial analysis of mortality in relation to BWQS scores ####
# Step 1: Create dataframes with the relevant information
deaths_by23May2020_by_zcta1 <- deaths_by23May2020_by_zcta %>%
left_join(., modzcta_to_zcta, by = c("MODIFIED_ZCTA" = "MODZCTA")) %>%
mutate(zcta = as.character(MODIFIED_ZCTA)) %>%
rename("deaths_count" = "COVID_DEATH_COUNT") %>%
dplyr::select(zcta, deaths_count, COVID_DEATH_RATE) %>%
distinct(zcta, .keep_all = T)
ZCTA_BWQS_COVID_shp1 <- ZCTA_ACS_COVID_shp %>%
left_join(.,ZCTA_BWQS_score, by = "zcta") %>%
left_join(., deaths_by23May2020_by_zcta1, by = "zcta") %>%
mutate(prop_65plus = age65_plus/total_pop1,
zcta = as.numeric(zcta))
# Figure 3C - Mortality
fig3c <- ZCTA_BWQS_COVID_shp1 %>%
ggplot() +
geom_sf(data = basemap_water, fill = "white", lwd = 0) +
geom_sf(aes(fill = COVID_DEATH_RATE), lwd = 0.2)+
scalebar(MODZCTA_NYC_shp1, dist = 5, dist_unit = "km",
transform = TRUE, model = "WGS84",
st.size = 2.8, height = 0.015, border.size = 0.5,
anchor = c(x = -73.71, y = 40.51)) +
labs(fill = "Mortality per 100,000") +
ggtitle("Cumulative COVID mortality by zip code (May 23, 2020)") +
scale_fill_gradientn(colours=brewer_pal("YlGnBu", type = "seq")(7)) +
coord_sf(crs = st_crs(ZCTA_BWQS_COVID_shp),
xlim = c(plot_bounds$xmin, plot_bounds$xmax),
ylim = c(plot_bounds$ymin, plot_bounds$ymax),
expand = FALSE) +
theme_bw(base_size = 6) +
theme_smallmaps
# Combine subfigures to make Figure 3: maps of BWQS, Tests, and Mortality by ZCTA
fig3a_2 = fig3a + labs(tag = "a") + theme(plot.tag.position = c(-0.018, 0.985), plot.tag = element_text(face = "bold", size = 15))
fig3b_2 = fig3b + labs(tag = "b") + theme(plot.tag.position = c(-0.028, 0.93), plot.tag = element_text(face = "bold", size = 15))
fig3c_2 = fig3c + labs(tag = "c") + theme(plot.tag.position = c(-0.028, 0.93), plot.tag = element_text(face = "bold", size = 15))
plotres = 300
agg_png(filename = file.path(fig.path, paste0("fig3_combined_", Sys.Date(), ".png")),
width = plotres*4.25, height = plotres*6, scaling = 2)
grid.arrange(arrangeGrob(fig3a_2, ncol = 1, nrow = 1),
arrangeGrob(fig3b_2, fig3c_2, ncol = 2, nrow = 1),
nrow = 2, ncol = 1, heights = c(1.9,1))
dev.off()## png
## 2
# Step 2: Run negative binomial model with spatial filtering
# Step 2a: Identify the neighbor relationships
spdat.sens <- as_Spatial(ZCTA_BWQS_COVID_shp1)
ny.nb6 <- knearneigh(sp::coordinates(spdat.sens), k=6)
ny.nb6 <- knn2nb(ny.nb6)
ny.nb6 <- make.sym.nb(ny.nb6)
ny.wt6 <- nb2listw(ny.nb6, style="W")
# Step 2b: Fit the model to identify the component of the data with substantial spatial autocorrelation
fit.nb.ny.sens<-glm.nb(deaths_count~offset(log(total_pop1))+BWQS_index, spdat.sens)
lm.morantest(fit.nb.ny.sens, listw = ny.wt6)##
## Global Moran I for regression residuals
##
## data:
## model: glm.nb(formula = deaths_count ~ offset(log(total_pop1)) +
## BWQS_index, data = spdat.sens, init.theta = 6.682813998, link = log)
## weights: ny.wt6
##
## Moran I statistic standard deviate = 3.1743, p-value = 0.000751
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.124763298 -0.001792051 0.001589528me.fit.sens <- spatialreg::ME(deaths_count~offset(log(total_pop1))+BWQS_index,
spdat.sens@data, family=negative.binomial(fit.nb.ny.sens$theta), listw = ny.wt6, verbose=T, alpha=.1, nsim = 999)## eV[,16], I: 0.05665305 ZI: NA, pr(ZI): 0.063
## eV[,23], I: 0.01207501 ZI: NA, pr(ZI): 0.342# Step 2c: Pull out these fits and visualize the autocorrelation
fits.sens <- data.frame(fitted(me.fit.sens))
spdat.sens$me16 <- fits.sens$vec16
spdat.sens$me23 <- fits.sens$vec23
spplot(spdat.sens, "me16", at=quantile(spdat.sens$me16, p=seq(0,1,length.out = 7)))
spplot(spdat.sens, "me23", at=quantile(spdat.sens$me23, p=seq(0,1,length.out = 7)))
# Step 2d: Include the fits in our regression model as an additional parameter
clean.nb.sens<-glm.nb(deaths_count~offset(log(total_pop1))+BWQS_index+fitted(me.fit.sens), spdat.sens@data)
tidy(clean.nb.sens) %>% mutate(estimate_exp = exp(estimate))## # A tibble: 4 x 6
## term estimate std.error statistic p.value estimate_exp
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) -7.29 0.0839 -86.9 0. 0.000679
## 2 BWQS_index 0.181 0.0151 12.0 3.73e-33 1.20
## 3 fitted(me.fit.sens)vec16 -0.925 0.383 -2.41 1.58e- 2 0.397
## 4 fitted(me.fit.sens)vec23 -1.78 0.381 -4.66 3.16e- 6 0.169as_tibble(confint(clean.nb.sens), rownames = "vars")%>% mutate_at(vars(2:3), .funs = list(~exp(.)))## Waiting for profiling to be done...## # A tibble: 4 x 3
## vars `2.5 %` `97.5 %`
## <chr> <dbl> <dbl>
## 1 (Intercept) 0.000575 0.000804
## 2 BWQS_index 1.16 1.23
## 3 fitted(me.fit.sens)vec16 0.188 0.836
## 4 fitted(me.fit.sens)vec23 0.0789 0.361lm.morantest(clean.nb.sens, resfun = residuals, listw=ny.wt6)##
## Global Moran I for regression residuals
##
## data:
## model: glm.nb(formula = deaths_count ~ offset(log(total_pop1)) +
## BWQS_index + fitted(me.fit.sens), data = spdat.sens@data, init.theta =
## 8.00423828, link = log)
## weights: ny.wt6
##
## Moran I statistic standard deviate = 1.643, p-value = 0.05019
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.063550312 -0.002549371 0.001618464spdat.sens$nb_residuals <- clean.nb.sens$residuals
spplot(spdat.sens, "nb_residuals", at=quantile(spdat.sens$nb_residuals, p=seq(0,1,length.out = 10)))
pal_sfig10 <- brewer_pal(type = "div", palette = "RdYlBu")(5)[2:5]
pal_sfig10[[2]] <- "#fafadf" # make yellow less saturated
sfig10 <- ggplot() +
geom_sf(data = basemap_water, fill = "white", lwd = 0) +
geom_sf(data = st_as_sf(spdat.sens), aes(fill = cut_width(nb_residuals, width = 1)), lwd = 0.2) +
scale_fill_discrete(type = pal_sfig10) +
coord_sf(crs = st_crs(MODZCTA_NYC_shp1),
xlim = c(plot_bounds$xmin, plot_bounds$xmax),
ylim = c(plot_bounds$ymin, plot_bounds$ymax),
expand = FALSE) +
labs(fill = "Residuals\n(Standard Deviations)") +
theme(panel.background = element_rect(fill = "#cccccc"),
panel.grid = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank(),
axis.title.y = element_blank(),
axis.title.x = element_blank(),
legend.title = element_text(face = "bold", size = 9),
legend.text = element_text(size = 8.5),
plot.margin = unit(c(0.1,0.1,0.1,0.1), "in")
)
if(export.figs) ggsave(sfig10, filename = file.path(fig.path, paste0("sfig10", "_", Sys.Date(),".png")),
dpi = 300, width = 5, height = 3.2)
# Sensitivity Analyses
#### Sensitivity Analyses ####
# half-Cauchy prior for overdispersion parameter
BWQS_NB_halfcauchy_stan_model <- here("code", "nb_bwqs_halfcauchy.stan")
m1_halfcauchy <- fit_BWQS_model(m1data$data_list, BWQS_NB_halfcauchy_stan_model)
# beta1 effect estimate using half-cauchy prior for overdispersion
exp(summary(m1_halfcauchy)$summary[3,c(1,4,8)])## mean 2.5% 97.5%
## 1.077273 1.066034 1.089665# main model (inverse gamma prior for overdispersion)
exp(summary(m1)$summary[3,c(1,4,8)])## mean 2.5% 97.5%
## 1.076975 1.064520 1.089281# no meaningful difference in effect estimate or credible interval!
# BWQS weights -- stability of the rankings
as_tibble(extract(m1, c("W[1]","W[2]", "W[3]", "W[4]", "W[5]", "W[6]", "W[7]", "W[8]", "W[9]", "W[10]"))) %>%
rownames_to_column("iteration") %>%
gather(weight, value, "W[1]":"W[10]") %>%
group_by(iteration) %>%
slice(which.max(value)) %>%
ungroup() %>%
group_by(weight) %>%
dplyr::summarise(times_largest_weight = n()) %>%
mutate(pct_largest_weight = round((times_largest_weight/1750)*100, 2))## # A tibble: 8 x 3
## weight times_largest_weight pct_largest_weight
## <chr> <int> <dbl>
## 1 W[1] 94 5.37
## 2 W[10] 457 26.1
## 3 W[2] 564 32.2
## 4 W[4] 4 0.23
## 5 W[5] 13 0.74
## 6 W[7] 504 28.8
## 7 W[8] 35 2
## 8 W[9] 79 4.51# using essential workers, testing ratio, and median income
glm_esstl_and_income <- glm.nb(pos_per_100000 ~ not_quarantined_jobs + medincome + testing_ratio, data = ZCTA_ACS_COVID)
broom::tidy(glm_esstl_and_income)## # A tibble: 4 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 6.34 0.0701 90.4 0.
## 2 not_quarantined_jobs 0.0185 0.00248 7.46 8.96e- 14
## 3 medincome -0.00000250 0.000000359 -6.96 3.36e- 12
## 4 testing_ratio 19.5 0.888 22.0 4.07e-107enframe(predict(glm_esstl_and_income)) %>% mutate(value = exp(value))## # A tibble: 177 x 2
## name value
## <chr> <dbl>
## 1 1 1293.
## 2 2 1301.
## 3 3 877.
## 4 4 870.
## 5 5 710.
## 6 6 787.
## 7 7 677.
## 8 8 1200.
## 9 9 987.
## 10 10 891.
## # … with 167 more rows# Just using proportion uninsured and median income
glm_insurance_and_income <- glm.nb(pos_per_100000 ~ not_insured + medincome, data = ZCTA_ACS_COVID)
broom::tidy(glm_insurance_and_income)## # A tibble: 3 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 8.14 0.120 68.0 0.
## 2 not_insured -0.00420 0.00857 -0.490 6.24e- 1
## 3 medincome -0.00000737 0.000000891 -8.26 1.41e-16# PCA of social variables
pca_socialvars <- prcomp(ZCTA_ACS_COVID %>% dplyr::select(all_of(SES_vars)), center = T, scale. = T)
PC1 <- enframe(pca_socialvars$rotation[,1]) %>% arrange(desc(value)) %>% rename("label" = "name")
df_for_PCA_analysis <- ZCTA_ACS_COVID %>% bind_cols(., PC1 = pca_socialvars$x[,1])
# Compare BWQS ranks to Principal Components ranks
BWQS_fits %>%
dplyr::select(mean, label) %>%
mutate(rank_bwqs = rank(mean)) %>%
left_join(., PC1, by = "label") %>%
mutate(rank_pca = rank(value)) %>%
dplyr::select(2, 3, 5, 1, 4)## # A tibble: 10 x 5
## label rank_bwqs rank_pca mean value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 not_insured 10 4 0.177 0.317
## 2 essentialworker_drove 9 2 0.166 0.201
## 3 avg_hhold_size 8 9 0.161 0.382
## 4 res_vol_popdensity 7 3 0.108 0.217
## 5 one_over_medincome 6 8 0.105 0.363
## 6 didnot_workhome_commute 5 7 0.0894 0.359
## 7 not_quarantined_jobs 4 10 0.0618 0.397
## 8 unemployed 3 6 0.0547 0.358
## 9 essentialworker_pubtrans 2 5 0.0498 0.335
## 10 one_over_grocers_per_1000 1 1 0.0267 0.0891glm_princomp <- glm.nb(pos_per_100000 ~ PC1 + ns(testing_ratio, df = 3), data = df_for_PCA_analysis)
summary(glm_princomp)##
## Call:
## glm.nb(formula = pos_per_100000 ~ PC1 + ns(testing_ratio, df = 3),
## data = df_for_PCA_analysis, init.theta = 94.1393225, link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5328 -0.5874 0.0735 0.5595 3.0695
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 6.635507 0.047025 141.10 <2e-16 ***
## PC1 0.071785 0.004996 14.37 <2e-16 ***
## ns(testing_ratio, df = 3)1 0.980565 0.042292 23.19 <2e-16 ***
## ns(testing_ratio, df = 3)2 2.033267 0.106250 19.14 <2e-16 ***
## ns(testing_ratio, df = 3)3 1.135913 0.053253 21.33 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(94.1393) family taken to be 1)
##
## Null deviance: 2797.59 on 176 degrees of freedom
## Residual deviance: 179.89 on 172 degrees of freedom
## AIC: 2381
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 94.1
## Std. Err.: 10.7
##
## 2 x log-likelihood: -2368.997exp(confint(glm.nb(pos_per_100000 ~ PC1 + ns(testing_ratio, df = 3), data = df_for_PCA_analysis)))## Waiting for profiling to be done...## 2.5 % 97.5 %
## (Intercept) 693.994592 836.504296
## PC1 1.063870 1.085062
## ns(testing_ratio, df = 3)1 2.456486 2.893115
## ns(testing_ratio, df = 3)2 6.178602 9.435865
## ns(testing_ratio, df = 3)3 2.805253 3.461891enframe(predict(glm_princomp)) %>% mutate(glm_princomp = exp(value)) %>% dplyr::select(glm_princomp)## # A tibble: 177 x 1
## glm_princomp
## <dbl>
## 1 1333.
## 2 1088.
## 3 727.
## 4 943.
## 5 741.
## 6 819.
## 7 891.
## 8 1065.
## 9 920.
## 10 809.
## # … with 167 more rowsenframe(predict(glm_esstl_and_income)) %>% mutate(glm_esstl_and_income = exp(value)) %>% dplyr::select(glm_esstl_and_income)## # A tibble: 177 x 1
## glm_esstl_and_income
## <dbl>
## 1 1293.
## 2 1301.
## 3 877.
## 4 870.
## 5 710.
## 6 787.
## 7 677.
## 8 1200.
## 9 987.
## 10 891.
## # … with 167 more rowsCompare_Metrics <- bind_cols(ZCTA_ACS_COVID,
glm_princomp = enframe(predict(glm_princomp)) %>% mutate(glm_princomp = exp(value)) %>% dplyr::select(glm_princomp),
glm_esstl_and_income = enframe(predict(glm_esstl_and_income)) %>% mutate(glm_esstl_and_income = exp(value)) %>% dplyr::select(glm_esstl_and_income),
BWQS_index = colMeans(extract(m1,"y_new")$y_new))
# comparison of RMSE and Kendall's tau metrics for three different modeling approaches
# Supplemental Table 4
bind_rows(Compare_Metrics %>%
summarise_at(vars(c("BWQS_index", "glm_princomp", "glm_esstl_and_income")),
~cor(pos_per_100000, .x, method = "kendall")) %>%
mutate(parameter = "kendalls"),
Compare_Metrics %>%
summarise_at(vars(c("BWQS_index", "glm_princomp", "glm_esstl_and_income")),
~sqrt(mean(pos_per_100000-.x)^2)) %>%
mutate(parameter = "RMSE"))## # A tibble: 2 x 4
## BWQS_index glm_princomp glm_esstl_and_income parameter
## <dbl> <dbl> <dbl> <chr>
## 1 0.869 0.850 0.823 kendalls
## 2 2.79 2.13 8.10 RMSESubway analyses
#### Subway analyses ####
# Different BWQS groupings for subway analysis (<.25, .25-.75, .75)
UHF_BWQS_COVID_3split_shp <- UHF_shp %>%
mutate(uhf = as.character(UHFCODE)) %>%
left_join(., UHF_BWQS, by = "uhf") %>%
mutate(Risk = factor(ifelse(uhf_weighted_bwqs < quantile(uhf_weighted_bwqs, probs = .25, na.rm = T), "Low",
if_else(uhf_weighted_bwqs > quantile(uhf_weighted_bwqs, probs = .75, na.rm = T), "High", "Mid")),
levels = c("High", "Mid", "Low")))
Subway_BWQS_3split_df <- Subway_ridership_by_UHF %>%
left_join(., st_set_geometry(UHF_BWQS_COVID_3split_shp, NULL), by = c("place" = "uhf")) %>%
filter(!is.na(place) & date >= "2020-01-29" & date <= "2020-04-30") %>%
group_by(place) %>%
arrange(date) %>%
mutate(time_index = time(date)) %>%
filter(!is.na(Risk) & date != "2020-02-17") %>% #removing Presidents' Day national holiday
anti_join(., service_changes_in_lowsubway_areas, by = c("date", "place")) #removing low outliers due to service changes in low subway density areas
fit_drm_risk_3split <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_3split_df, separate = F)
anova(fit_drm_all, fit_drm_risk_3split)##
## 1st model
## fct: W2.4()
## pmodels: 1 (for all parameters)
## 2nd model
## fct: W2.4()
## pmodels: Risk (for all parameters)## ANOVA table
##
## ModelDf RSS Df F value p value
## 1st model 3291 21.302
## 2nd model 3283 15.432 8 156.1 0.0summary(fit_drm_risk_3split)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:Low -12.0668117 0.4809539 -25.089 < 2.2e-16 ***
## b:Mid -11.0653059 0.3525339 -31.388 < 2.2e-16 ***
## b:High -10.8697531 0.3999555 -27.177 < 2.2e-16 ***
## c:Low 0.0728584 0.0038983 18.690 < 2.2e-16 ***
## c:Mid 0.1294325 0.0033939 38.137 < 2.2e-16 ***
## c:High 0.1751494 0.0040700 43.034 < 2.2e-16 ***
## d:Low 0.9805997 0.0035709 274.612 < 2.2e-16 ***
## d:Mid 0.9598153 0.0028993 331.053 < 2.2e-16 ***
## [ reached getOption("max.print") -- omitted 4 rows ]
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.06856061 (3283 degrees of freedom)compParm(fit_drm_risk_3split, "b", operator = "-")##
## Comparison of parameter 'b'
##
## Estimate Std. Error t-value p-value
## Low-Mid -1.00151 0.59632 -1.6795 0.09315 .
## Low-High -1.19706 0.62552 -1.9137 0.05575 .
## Mid-High -0.19555 0.53315 -0.3668 0.71380
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1compParm(fit_drm_risk_3split, "c", operator = "-")##
## Comparison of parameter 'c'
##
## Estimate Std. Error t-value p-value
## Low-Mid -0.0565742 0.0051687 -10.9456 < 2.2e-16 ***
## Low-High -0.1022910 0.0056357 -18.1504 < 2.2e-16 ***
## Mid-High -0.0457168 0.0052994 -8.6268 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(fit_drm_risk_3split)
fit_drm_risk_3split_high <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_3split_df, subset = Risk=="High")
fit_drm_risk_3split_mid <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_3split_df, subset = Risk=="Mid")
fit_drm_risk_3split_low <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_3split_df, subset = Risk=="Low")
summary(fit_drm_risk_3split_high)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -10.8693720 0.3278018 -33.158 < 2.2e-16 ***
## c:(Intercept) 0.1751430 0.0033359 52.502 < 2.2e-16 ***
## d:(Intercept) 0.9758825 0.0027014 361.256 < 2.2e-16 ***
## e:(Intercept) 46.8619040 0.1050027 446.292 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.05619323 (1008 degrees of freedom)summary(fit_drm_risk_3split_mid)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -11.0658766 0.3964936 -27.909 < 2.2e-16 ***
## c:(Intercept) 0.1294362 0.0038167 33.913 < 2.2e-16 ***
## d:(Intercept) 0.9598145 0.0032605 294.376 < 2.2e-16 ***
## e:(Intercept) 45.9159522 0.1187447 386.678 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.07710264 (1359 degrees of freedom)summary(fit_drm_risk_3split_low)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -12.0672287 0.4733213 -25.495 < 2.2e-16 ***
## c:(Intercept) 0.0728676 0.0038363 18.994 < 2.2e-16 ***
## d:(Intercept) 0.9805959 0.0035141 279.045 < 2.2e-16 ***
## e:(Intercept) 44.2637382 0.1136008 389.643 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.06747137 (916 degrees of freedom)confint(fit_drm_risk_3split_high)## 2.5 % 97.5 %
## b:(Intercept) -11.5126240 -10.2261199
## c:(Intercept) 0.1685968 0.1816891
## d:(Intercept) 0.9705815 0.9811834
## e:(Intercept) 46.6558551 47.0679529confint(fit_drm_risk_3split_mid)## 2.5 % 97.5 %
## b:(Intercept) -11.8436826 -10.2880707
## c:(Intercept) 0.1219489 0.1369234
## d:(Intercept) 0.9534184 0.9662107
## e:(Intercept) 45.6830095 46.1488949confint(fit_drm_risk_3split_low)## 2.5 % 97.5 %
## b:(Intercept) -12.99614880 -11.13830857
## c:(Intercept) 0.06533873 0.08039652
## d:(Intercept) 0.97369930 0.98749260
## e:(Intercept) 44.04079005 44.48668636# b is not actually a slope - but a scaling factor that needs to be transformed into the slope
slopes_3split <- as_tibble(coef(fit_drm_risk_3split_high), rownames = "vars") %>%
bind_cols(BWQS_risk = "high") %>%
bind_rows(as_tibble(coef(fit_drm_risk_3split_mid), rownames = "vars") %>% bind_cols(BWQS_risk = "mid")) %>%
bind_rows(as_tibble(coef(fit_drm_risk_3split_low), rownames = "vars") %>% bind_cols(BWQS_risk = "low")) %>%
mutate(vars = str_remove(vars, fixed(":(Intercept)"))) %>%
spread(., vars, value) %>%
mutate(slope = -b*(-d/(4*e)))
as_tibble(confint(fit_drm_risk_3split_high), rownames = "vars") %>%
bind_cols(BWQS_risk = "high") %>%
bind_rows(as_tibble(confint(fit_drm_risk_3split_mid), rownames = "vars") %>% bind_cols(BWQS_risk = "mid")) %>%
bind_rows(as_tibble(confint(fit_drm_risk_3split_low), rownames = "vars") %>% bind_cols(BWQS_risk = "low")) %>%
mutate(vars = str_remove(vars, fixed(":(Intercept)"))) %>%
filter(vars == "b") %>%
right_join(., slopes_3split, by = "BWQS_risk") %>%
mutate(slope_lower_ci = -`2.5 %`*(-d/(4*e)),
slope_higher_ci = -`97.5 %`*(-d/(4*e))) %>%
dplyr::select(BWQS_risk, slope, slope_lower_ci, slope_higher_ci)## # A tibble: 3 x 4
## BWQS_risk slope slope_lower_ci slope_higher_ci
## <chr> <dbl> <dbl> <dbl>
## 1 high -0.0566 -0.0599 -0.0532
## 2 mid -0.0578 -0.0619 -0.0538
## 3 low -0.0668 -0.0720 -0.0617fit_drm_predictions_3split <- as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_high, interval = "confidence"), warning = handler) %>%
bind_cols(., Risk = "High")) %>%
bind_rows(., as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_mid, interval = "confidence"),
warning = handler ) %>% bind_cols(., Risk = "Mid"))) %>%
bind_rows(., as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_low, interval = "confidence"),
warning = handler ) %>% bind_cols(., Risk = "Low")))
Subway_BWQS_3split_df1 <- bind_cols(Subway_BWQS_3split_df %>% arrange(Risk, date),
fit_drm_predictions_3split %>% dplyr::select(-Risk))
Subway_BWQS_3split_df2 <- Subway_BWQS_3split_df1 %>%
filter(date>"2020-02-16") %>% # subsetting for visualization
mutate(Risk = factor(if_else(Risk == "High", "High (≥ 75%ile)",
if_else(Risk=="Mid", "Mid (IQR)", "Low (≤ 25%ile)")),
levels = c("High (≥ 75%ile)", "Mid (IQR)", "Low (≤ 25%ile)")))
# Supplementary Figure 7
sfig7 <- ggplot() +
geom_jitter(data = Subway_BWQS_3split_df2, aes(x = date, y = usage.median.ratio, color = as.factor(Risk)), alpha = .5, position = position_jitter(height = 0, width = 0.4))+
geom_ribbon(data = subset(Subway_BWQS_3split_df2, Risk == "High (≥ 75%ile)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_ribbon(data = subset(Subway_BWQS_3split_df2, Risk == "Mid (IQR)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_ribbon(data = subset(Subway_BWQS_3split_df2, Risk == "Low (≤ 25%ile)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_line(data = Subway_BWQS_3split_df2, aes(x = date, y = Prediction, color = as.factor(Risk))) +
scale_x_date("Date", date_minor_breaks = "1 week") +
scale_y_continuous("Relative Subway Ridership (%)", labels = scales::percent) +
geom_vline(xintercept = date("2020-03-22"),color = "grey30", lty = 2) +
theme_bw(base_size = 16) +
labs(colour="COVID-19 inequity index") +
theme(legend.title = element_text(face = "bold", size = 12), legend.position = c(0.8, 0.7))
if(export.figs) ggsave(sfig7 , filename = file.path(fig.path, paste0("sfig7 ", "_", Sys.Date(),".png")),
dpi = 300, width = 8, height = 6)
# ZCTA level BWQS for subway
Subway_ridership_by_ZCTA <- relative.subway.usage(2020, by = "zcta")
# There are clearly some higher outliers that need to be removed at the zcta level
# Remove the zctas associated with the uhfs below
high_april_bronx_subway <- tibble(date = as.Date(c("2020-04-04", "2020-04-05", "2020-04-11", "2020-04-12", "2020-04-18", "2020-04-19", "2020-04-25", "2020-04-26")),
place = c("10469", "10469","10469","10469","10469","10469","10469","10469"))
Subway_BWQS_ZCTA_df <- Subway_ridership_by_ZCTA %>%
left_join(., ZCTA_BWQS_score, by = c("place" = "zcta")) %>%
filter(!is.na(place) & date >= "2020-01-29" & date <= "2020-04-30") %>%
group_by(place) %>%
arrange(date) %>%
mutate(time_index = time(date)) %>%
ungroup() %>%
filter(!is.na(BWQS_index) & date != "2020-02-17") %>% #removing Presidents' Day national holiday
anti_join(., high_april_bronx_subway, by = c("date", "place")) %>% #removing low outliers due to service changes in low subway density areas
mutate(Risk = factor(if_else(BWQS_index < quantile(BWQS_index, probs = .25), "Low",
if_else(BWQS_index > quantile(BWQS_index, probs = .75), "High", "Mid")),
levels = c("High", "Mid", "Low")))
fit_drm_risk_3split_zcta_high <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_ZCTA_df, subset = Risk=="High")
fit_drm_risk_3split_zcta_mid <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_ZCTA_df, subset = Risk=="Mid")
fit_drm_risk_3split_zcta_low <- drm(usage.median.ratio ~ time_index, curveid = Risk, fct = W2.4(), data = Subway_BWQS_ZCTA_df, subset = Risk=="Low")
summary(fit_drm_risk_3split_zcta_high)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -10.9137233 0.3404024 -32.061 < 2.2e-16 ***
## c:(Intercept) 0.1686448 0.0034263 49.221 < 2.2e-16 ***
## d:(Intercept) 0.9566853 0.0027549 347.266 < 2.2e-16 ***
## e:(Intercept) 46.9474325 0.1081606 434.053 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.09486494 (2748 degrees of freedom)summary(fit_drm_risk_3split_zcta_mid)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -11.3113739 0.3758877 -30.092 < 2.2e-16 ***
## c:(Intercept) 0.1340174 0.0035621 37.623 < 2.2e-16 ***
## d:(Intercept) 0.9817241 0.0030165 325.456 < 2.2e-16 ***
## e:(Intercept) 45.9104803 0.1072824 427.940 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.1474759 (5700 degrees of freedom)summary(fit_drm_risk_3split_zcta_low)##
## Model fitted: Weibull (type 2) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -11.8046717 0.3455007 -34.167 < 2.2e-16 ***
## c:(Intercept) 0.0619258 0.0028874 21.447 < 2.2e-16 ***
## d:(Intercept) 0.9886886 0.0026770 369.327 < 2.2e-16 ***
## e:(Intercept) 43.9322297 0.0871923 503.854 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 0.08812689 (2756 degrees of freedom)# b is not actually a slope - but a scaling factor that needs to be transformed into the slope
slopes_zcta <- as_tibble(coef(fit_drm_risk_3split_zcta_high), rownames = "vars") %>%
bind_cols(BWQS_risk = "high") %>%
bind_rows(as_tibble(coef(fit_drm_risk_3split_zcta_mid), rownames = "vars") %>% bind_cols(BWQS_risk = "mid")) %>%
bind_rows(as_tibble(coef(fit_drm_risk_3split_zcta_low), rownames = "vars") %>% bind_cols(BWQS_risk = "low")) %>%
mutate(vars = str_remove(vars, fixed(":(Intercept)"))) %>%
spread(., vars, value) %>%
mutate(slope = -b*(-d/(4*e)))
as_tibble(confint(fit_drm_risk_3split_zcta_high), rownames = "vars") %>%
bind_cols(BWQS_risk = "high") %>%
bind_rows(as_tibble(confint(fit_drm_risk_3split_zcta_mid), rownames = "vars") %>% bind_cols(BWQS_risk = "mid")) %>%
bind_rows(as_tibble(confint(fit_drm_risk_3split_zcta_low), rownames = "vars") %>% bind_cols(BWQS_risk = "low")) %>%
mutate(vars = str_remove(vars, fixed(":(Intercept)"))) %>%
filter(vars == "b") %>%
right_join(., slopes_zcta, by = "BWQS_risk") %>%
mutate(slope_lower_ci = -`2.5 %`*(-d/(4*e)),
slope_higher_ci = -`97.5 %`*(-d/(4*e))) %>%
dplyr::select(BWQS_risk, slope, slope_lower_ci, slope_higher_ci)## # A tibble: 3 x 4
## BWQS_risk slope slope_lower_ci slope_higher_ci
## <chr> <dbl> <dbl> <dbl>
## 1 high -0.0556 -0.0590 -0.0522
## 2 mid -0.0605 -0.0644 -0.0565
## 3 low -0.0664 -0.0702 -0.0626fit_drm_predictions_zcta <- as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_zcta_high, interval = "confidence"), warning = handler) %>%
bind_cols(., Risk = "High")) %>%
bind_rows(., as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_zcta_mid, interval = "confidence"),
warning = handler ) %>% bind_cols(., Risk = "Mid"))) %>%
bind_rows(., as_tibble(withCallingHandlers(predict(fit_drm_risk_3split_zcta_low, interval = "confidence"),
warning = handler ) %>% bind_cols(., Risk = "Low"))) # %>% mutate(time_index = row_number())
Subway_BWQS_ZCTA_df1 <- bind_cols(Subway_BWQS_ZCTA_df %>% arrange(Risk, date),
fit_drm_predictions_zcta %>% dplyr::select(-Risk))
Subway_BWQS_ZCTA_df2 <- Subway_BWQS_ZCTA_df1 %>%
filter(date>"2020-02-16") %>% # subsetting for visualization
mutate(Risk = factor(if_else(Risk == "High", "High (≥ 75%ile)",
if_else(Risk=="Mid", "Mid (IQR)", "Low (≤ 25%ile)")),
levels = c("High (≥ 75%ile)", "Mid (IQR)", "Low (≤ 25%ile)")))
# Supplementary Figure 8
sfig8 <- ggplot() +
geom_jitter(data = Subway_BWQS_ZCTA_df2, aes(x = date, y = usage.median.ratio, color = Risk), size = 0.3, alpha = 0.3, position = position_jitter(height = 0, width = 0.42))+
geom_ribbon(data = subset(Subway_BWQS_ZCTA_df2, Risk == "High (≥ 75%ile)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_ribbon(data = subset(Subway_BWQS_ZCTA_df2, Risk == "Mid (IQR)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_ribbon(data = subset(Subway_BWQS_ZCTA_df2, Risk == "Low (≤ 25%ile)"), aes(x = date, ymin = Lower, ymax = Upper), fill = "grey50") +
geom_line(data = Subway_BWQS_ZCTA_df2, aes(x = date, y = Prediction, color = Risk)) +
scale_x_date("Date", date_minor_breaks = "1 week") +
scale_y_continuous("Relative Subway Ridership (%)", labels = scales::percent) +
geom_vline(xintercept = date("2020-03-22"),color = "grey30", lty = 2) +
theme_bw(base_size = 16) +
labs(colour="COVID-19 inequity index") +
theme(legend.title = element_text(face = "bold", size = 12), legend.position = c(0.8, 0.7))
if(export.figs) ggsave(sfig8 , filename = file.path(fig.path, paste0("sfig8", "_", Sys.Date(),".png")),
dpi = 300, width = 8, height = 6)
# Supplementary Figure 9: compare MTA turnstile data to Google mobility reportssource(here("code/mta_vs_google.R"), echo = FALSE)
mobplot
Appendix
#### Appendix ####
t_final = Sys.time()Runtimes
print(t_start)## [1] "2020-11-24 19:37:44 EST"print(t_turnstile_1)## [1] "2020-11-24 19:37:49 EST"print(t_turnstile_2)## [1] "2020-11-24 19:52:05 EST"print(t_census)## [1] "2020-11-24 19:57:07 EST"print(t_dataprep)## [1] "2020-11-24 20:05:43 EST"print(t_m1)## [1] "2020-11-24 20:13:43 EST"print(t_postm1)## [1] "2020-11-24 20:13:52 EST"print(t_m2)## [1] "2020-11-24 20:20:06 EST"print(t_m3)## [1] "2020-11-24 20:24:50 EST"print(t_final)## [1] "2020-11-24 20:32:08 EST"# Total runtime:
print(t_final - t_start)## Time difference of 54.40026 mins# Time part 1: Script start through loading of MTA Turnstile package:
print(t_turnstile_1 - t_start)## Time difference of 5.607934 secs# Time part 2: Downloading (if this is the first run) and processing of subway turnstile data by neighborhood:
print(t_turnstile_2 - t_turnstile_1)## Time difference of 14.25564 mins# Time part 3: Downloading of all other data and processing of Census data:
print(t_census - t_turnstile_2)## Time difference of 5.045368 mins# Time part 4: All other data processing, modeling, and figure generation:
print(t_final - t_census)## Time difference of 35.00579 minsSession Info
sessioninfo::session_info()## ─ Session info ───────────────────────────────────────────────────────────────────────────────────
## setting value
## version R version 4.0.2 (2020-06-22)
## os Ubuntu 18.04.4 LTS
## system x86_64, linux-gnu
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/New_York
## date 2020-11-24
##
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────
## package * version date lib source
## abind * 1.4-5 2016-07-21 [2] CRAN (R 4.0.2)
## acepack 1.4.1 2016-10-29 [2] CRAN (R 4.0.2)
## askpass 1.1 2019-01-13 [2] CRAN (R 4.0.2)
## assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.0.2)
## backports 1.1.8 2020-06-17 [2] CRAN (R 4.0.2)
## base64enc 0.1-3 2015-07-28 [2] CRAN (R 4.0.2)
## BBmisc 1.11 2017-03-10 [2] CRAN (R 4.0.2)
## bit 1.1-15.2 2020-02-10 [2] CRAN (R 4.0.2)
## bit64 0.9-7 2017-05-08 [2] CRAN (R 4.0.2)
## bitops 1.0-6 2013-08-17 [2] CRAN (R 4.0.2)
## blob 1.2.1 2020-01-20 [2] CRAN (R 4.0.2)
## boot 1.3-25 2020-04-26 [3] CRAN (R 4.0.0)
## broom * 0.5.6 2020-04-20 [2] CRAN (R 4.0.2)
## callr 3.4.3 2020-03-28 [2] CRAN (R 4.0.2)
## car 3.0-8 2020-05-21 [2] CRAN (R 4.0.2)
## carData 3.0-4 2020-05-22 [2] CRAN (R 4.0.2)
## cellranger 1.1.0 2016-07-27 [2] CRAN (R 4.0.2)
## checkmate 2.0.0 2020-02-06 [2] CRAN (R 4.0.2)
## class 7.3-17 2020-04-26 [3] CRAN (R 4.0.0)
## classInt 0.4-3 2020-04-07 [2] CRAN (R 4.0.2)
## cli 2.0.2 2020-02-28 [2] CRAN (R 4.0.2)
## cluster 2.1.0 2019-06-19 [3] CRAN (R 4.0.0)
## coda 0.19-3 2019-07-05 [2] CRAN (R 4.0.2)
## codetools 0.2-16 2018-12-24 [3] CRAN (R 4.0.0)
## colorspace 1.4-1 2019-03-18 [2] CRAN (R 4.0.2)
## corrplot * 0.84 2017-10-16 [2] CRAN (R 4.0.2)
## crayon 1.3.4 2017-09-16 [2] CRAN (R 4.0.2)
## curl 4.3 2019-12-02 [2] CRAN (R 4.0.2)
## data.table * 1.13.0 2020-07-24 [1] CRAN (R 4.0.2)
## DBI 1.1.0 2019-12-15 [2] CRAN (R 4.0.2)
## dbplyr 1.4.4 2020-05-27 [2] CRAN (R 4.0.2)
## deldir 0.1-25 2020-02-03 [2] CRAN (R 4.0.2)
## DHARMa * 0.3.3.0 2020-09-08 [2] CRAN (R 4.0.2)
## digest 0.6.25 2020-02-23 [2] CRAN (R 4.0.2)
## dplyr * 1.0.0 2020-05-29 [2] CRAN (R 4.0.2)
## drc * 3.0-1 2016-08-30 [1] CRAN (R 4.0.2)
## e1071 1.7-3 2019-11-26 [2] CRAN (R 4.0.2)
## egg * 0.4.5 2019-07-13 [2] CRAN (R 4.0.2)
## ellipsis 0.3.1 2020-05-15 [2] CRAN (R 4.0.2)
## evaluate 0.14 2019-05-28 [2] CRAN (R 4.0.2)
## expm 0.999-4 2019-03-21 [2] CRAN (R 4.0.2)
## fansi 0.4.1 2020-01-08 [2] CRAN (R 4.0.2)
## farver 2.0.3 2020-01-16 [2] CRAN (R 4.0.2)
## fastmap 1.0.1 2019-10-08 [2] CRAN (R 4.0.2)
## fastmatch 1.1-0 2017-01-28 [2] CRAN (R 4.0.2)
## forcats * 0.5.0 2020-03-01 [2] CRAN (R 4.0.2)
## foreach 1.5.0 2020-03-30 [2] CRAN (R 4.0.2)
## foreign 0.8-79 2020-04-26 [3] CRAN (R 4.0.0)
## Formula 1.2-3 2018-05-03 [2] CRAN (R 4.0.2)
## fs 1.4.2 2020-06-30 [2] CRAN (R 4.0.2)
## fst * 0.9.2 2020-04-01 [2] CRAN (R 4.0.2)
## gap 1.2.2 2020-02-02 [2] CRAN (R 4.0.2)
## gdata 2.18.0 2017-06-06 [2] CRAN (R 4.0.2)
## generics 0.0.2 2018-11-29 [2] CRAN (R 4.0.2)
## ggExtra * 0.9 2019-08-27 [2] CRAN (R 4.0.2)
## ggmap 3.0.0 2019-02-05 [2] CRAN (R 4.0.2)
## ggplot2 * 3.3.2 2020-06-19 [2] CRAN (R 4.0.2)
## ggpubr * 0.4.0 2020-06-27 [1] CRAN (R 4.0.2)
## ggridges * 0.5.2 2020-01-12 [2] CRAN (R 4.0.2)
## ggsignif 0.6.0 2019-08-08 [1] CRAN (R 4.0.2)
## ggsn * 0.5.0 2019-02-18 [2] CRAN (R 4.0.2)
## glue 1.4.2 2020-08-27 [1] CRAN (R 4.0.2)
## gmodels 2.18.1 2018-06-25 [2] CRAN (R 4.0.2)
## gridExtra * 2.3 2017-09-09 [2] CRAN (R 4.0.2)
## gtable 0.3.0 2019-03-25 [2] CRAN (R 4.0.2)
## gtools 3.8.2 2020-03-31 [2] CRAN (R 4.0.2)
## haven 2.3.1 2020-06-01 [2] CRAN (R 4.0.2)
## here * 0.1 2017-05-28 [2] CRAN (R 4.0.2)
## highr 0.8 2019-03-20 [2] CRAN (R 4.0.2)
## Hmisc 4.4-0 2020-03-23 [2] CRAN (R 4.0.2)
## hms 0.5.3 2020-01-08 [2] CRAN (R 4.0.2)
## htmlTable 2.0.0 2020-06-21 [2] CRAN (R 4.0.2)
## htmltools 0.5.0 2020-06-16 [2] CRAN (R 4.0.2)
## htmlwidgets 1.5.1 2019-10-08 [2] CRAN (R 4.0.2)
## httpuv 1.5.4 2020-06-06 [2] CRAN (R 4.0.2)
## httr * 1.4.1 2019-08-05 [2] CRAN (R 4.0.2)
## inline 0.3.15 2018-05-18 [2] CRAN (R 4.0.2)
## iterators 1.0.12 2019-07-26 [2] CRAN (R 4.0.2)
## jpeg 0.1-8.1 2019-10-24 [2] CRAN (R 4.0.2)
## jsonlite * 1.7.0 2020-06-25 [2] CRAN (R 4.0.2)
## Just.universal * 0.0.0.9000 2020-11-24 [1] Github (justlab/Just_universal@4e6147a)
## kableExtra * 1.3.1 2020-10-22 [1] CRAN (R 4.0.2)
## KernSmooth 2.23-17 2020-04-26 [3] CRAN (R 4.0.0)
## knitr 1.29 2020-06-23 [2] CRAN (R 4.0.2)
## labeling 0.3 2014-08-23 [2] CRAN (R 4.0.2)
## later 1.1.0.1 2020-06-05 [2] CRAN (R 4.0.2)
## lattice 0.20-41 2020-04-02 [3] CRAN (R 4.0.0)
## latticeExtra 0.6-29 2019-12-19 [2] CRAN (R 4.0.2)
## LearnBayes 2.15.1 2018-03-18 [2] CRAN (R 4.0.2)
## lifecycle 0.2.0 2020-03-06 [2] CRAN (R 4.0.2)
## lme4 1.1-23 2020-04-07 [2] CRAN (R 4.0.2)
## loo 2.2.0 2019-12-19 [2] CRAN (R 4.0.2)
## lubridate * 1.7.9 2020-06-08 [2] CRAN (R 4.0.2)
## lwgeom * 0.2-5 2020-06-12 [2] CRAN (R 4.0.2)
## magic * 1.5-9 2018-09-17 [2] CRAN (R 4.0.2)
## magrittr 1.5 2014-11-22 [2] CRAN (R 4.0.2)
## maptools 1.0-1 2020-05-14 [2] CRAN (R 4.0.2)
## MASS * 7.3-51.6 2020-04-26 [3] CRAN (R 4.0.0)
## Matrix * 1.2-18 2019-11-27 [3] CRAN (R 4.0.0)
## matrixStats * 0.56.0 2020-03-13 [2] CRAN (R 4.0.2)
## memoise 1.1.0 2017-04-21 [2] CRAN (R 4.0.2)
## mgcv 1.8-31 2019-11-09 [3] CRAN (R 4.0.0)
## mime 0.9 2020-02-04 [2] CRAN (R 4.0.2)
## miniUI 0.1.1.1 2018-05-18 [2] CRAN (R 4.0.2)
## minqa 1.2.4 2014-10-09 [2] CRAN (R 4.0.2)
## modelr 0.1.8 2020-05-19 [2] CRAN (R 4.0.2)
## MTA.turnstile * 0.0.0.9001 2020-11-24 [1] Github (justlab/MTA_turnstile@544eba3)
## multcomp 1.4-13 2020-04-08 [2] CRAN (R 4.0.2)
## munsell 0.5.0 2018-06-12 [2] CRAN (R 4.0.2)
## mvtnorm 1.1-1 2020-06-09 [2] CRAN (R 4.0.2)
## nlme 3.1-147 2020-04-13 [3] CRAN (R 4.0.0)
## nloptr 1.2.2.1 2020-03-11 [2] CRAN (R 4.0.2)
## nnet 7.3-14 2020-04-26 [3] CRAN (R 4.0.0)
## openxlsx 4.1.5 2020-05-06 [2] CRAN (R 4.0.2)
## ParamHelpers 1.14 2020-03-24 [2] CRAN (R 4.0.2)
## pbapply 1.4-3 2020-08-18 [1] CRAN (R 4.0.2)
## pdftools * 2.3.1 2020-05-22 [1] CRAN (R 4.0.2)
## pillar 1.4.4 2020-05-05 [2] CRAN (R 4.0.2)
## pkgbuild 1.0.8 2020-05-07 [2] CRAN (R 4.0.2)
## pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.0.2)
## plotrix 3.7-8 2020-04-16 [2] CRAN (R 4.0.2)
## plyr 1.8.6 2020-03-03 [2] CRAN (R 4.0.2)
## png 0.1-7 2013-12-03 [2] CRAN (R 4.0.2)
## prettyunits 1.1.1 2020-01-24 [2] CRAN (R 4.0.2)
## processx 3.4.2 2020-02-09 [2] CRAN (R 4.0.2)
## promises 1.1.1 2020-06-09 [2] CRAN (R 4.0.2)
## ps 1.3.3 2020-05-08 [2] CRAN (R 4.0.2)
## purrr * 0.3.4 2020-04-17 [2] CRAN (R 4.0.2)
## qpdf 1.1 2019-03-07 [1] CRAN (R 4.0.2)
## qs * 0.22.1 2020-06-10 [2] CRAN (R 4.0.2)
## R6 2.4.1 2019-11-12 [2] CRAN (R 4.0.2)
## ragg * 0.4.0 2020-10-05 [1] CRAN (R 4.0.2)
## RANN 2.6.1 2019-01-08 [2] CRAN (R 4.0.2)
## RApiSerialize 0.1.0 2014-04-19 [2] CRAN (R 4.0.2)
## rappdirs 0.3.1 2016-03-28 [2] CRAN (R 4.0.2)
## raster 3.3-7 2020-06-27 [2] CRAN (R 4.0.2)
## RColorBrewer 1.1-2 2014-12-07 [2] CRAN (R 4.0.2)
## Rcpp 1.0.4.6 2020-04-09 [2] CRAN (R 4.0.2)
## RcppParallel 5.0.2 2020-06-24 [2] CRAN (R 4.0.2)
## readr * 1.3.1 2018-12-21 [2] CRAN (R 4.0.2)
## readxl * 1.3.1 2019-03-13 [2] CRAN (R 4.0.2)
## reprex 0.3.0 2019-05-16 [2] CRAN (R 4.0.2)
## rgdal 1.5-12 2020-06-26 [2] CRAN (R 4.0.2)
## RgoogleMaps 1.4.5.3 2020-02-12 [2] CRAN (R 4.0.2)
## rio 0.5.16 2018-11-26 [2] CRAN (R 4.0.2)
## rjson 0.2.20 2018-06-08 [2] CRAN (R 4.0.2)
## rlang 0.4.6 2020-05-02 [2] CRAN (R 4.0.2)
## rmarkdown 2.3 2020-06-18 [2] CRAN (R 4.0.2)
## rpart 4.1-15 2019-04-12 [3] CRAN (R 4.0.0)
## rprojroot 1.3-2 2018-01-03 [2] CRAN (R 4.0.2)
## RSQLite 2.2.1 2020-09-30 [2] CRAN (R 4.0.2)
## rstan * 2.19.3 2020-02-11 [2] CRAN (R 4.0.2)
## rstatix 0.6.0 2020-06-18 [1] CRAN (R 4.0.2)
## rstudioapi 0.11 2020-02-07 [2] CRAN (R 4.0.2)
## rvest 0.3.5 2019-11-08 [2] CRAN (R 4.0.2)
## sandwich 2.5-1 2019-04-06 [2] CRAN (R 4.0.2)
## scales * 1.1.1 2020-05-11 [2] CRAN (R 4.0.2)
## sessioninfo 1.1.1 2018-11-05 [2] CRAN (R 4.0.2)
## sf * 0.9-6 2020-09-13 [2] CRAN (R 4.0.2)
## shiny 1.5.0 2020-06-23 [2] CRAN (R 4.0.2)
## sp * 1.4-2 2020-05-20 [2] CRAN (R 4.0.2)
## spatialreg * 1.1-5 2019-12-01 [1] CRAN (R 4.0.2)
## spData * 0.3.5 2020-04-06 [2] CRAN (R 4.0.2)
## spdep * 1.1-5 2020-06-29 [2] CRAN (R 4.0.2)
## StanHeaders * 2.21.0-6 2020-08-16 [1] CRAN (R 4.0.2)
## statmod 1.4.34 2020-02-17 [2] CRAN (R 4.0.2)
## stringfish 0.12.1 2020-06-05 [2] CRAN (R 4.0.2)
## stringi 1.5.3 2020-09-09 [1] CRAN (R 4.0.2)
## stringr * 1.4.0 2019-02-10 [2] CRAN (R 4.0.2)
## survival 3.1-12 2020-04-10 [3] CRAN (R 4.0.0)
## systemfonts 0.3.2 2020-09-29 [1] CRAN (R 4.0.2)
## textshaping 0.1.2 2020-10-08 [1] CRAN (R 4.0.2)
## TH.data 1.0-10 2019-01-21 [2] CRAN (R 4.0.2)
## tibble * 3.0.1 2020-04-20 [2] CRAN (R 4.0.2)
## tidycensus * 0.10.2 2020-09-28 [1] CRAN (R 4.0.2)
## tidyr * 1.1.0 2020-05-20 [2] CRAN (R 4.0.2)
## tidyselect 1.1.0 2020-05-11 [2] CRAN (R 4.0.2)
## tidyverse * 1.3.0 2019-11-21 [2] CRAN (R 4.0.2)
## tigris 0.9.4 2020-04-01 [2] CRAN (R 4.0.2)
## units 0.6-7 2020-06-13 [2] CRAN (R 4.0.2)
## utf8 1.1.4 2018-05-24 [2] CRAN (R 4.0.2)
## uuid 0.1-4 2020-02-26 [2] CRAN (R 4.0.2)
## vctrs 0.3.1 2020-06-05 [2] CRAN (R 4.0.2)
## viridisLite 0.3.0 2018-02-01 [2] CRAN (R 4.0.2)
## webshot 0.5.2 2019-11-22 [2] CRAN (R 4.0.2)
## withr 2.2.0 2020-04-20 [2] CRAN (R 4.0.2)
## xfun 0.15 2020-06-21 [2] CRAN (R 4.0.2)
## xgboost 1.1.1.1 2020-06-14 [2] CRAN (R 4.0.2)
## xml2 1.3.2 2020-04-23 [2] CRAN (R 4.0.2)
## xtable 1.8-4 2019-04-21 [2] CRAN (R 4.0.2)
## yaml 2.2.1 2020-02-01 [2] CRAN (R 4.0.2)
## zip * 2.0.4 2019-09-01 [2] CRAN (R 4.0.2)
## zoo 1.8-8 2020-05-02 [2] CRAN (R 4.0.2)
##
## [1] /home/rushj03/R/x86_64-pc-linux-gnu-library/4.0
## [2] /usr/lib/R/site-library
## [3] /usr/lib/R/library